{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "This notebook gives suggests how to solve the problem of non-linear compressible flow using the automatic differentiation library included in PorePy. \n",
    "\n",
    "The Ad functionality in PorePy is currently undergoing a substantial expansion, see the tutorial AdFramework for indications of possibilities beyond what is covered herein."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model\n",
    "As an example, we will set up a non-linear problem for compressible flow. As usuall, we assume Darcy's law is valid:\n",
    "$$\n",
    "\\vec u = -\\mathcal K \\nabla p,\n",
    "$$\n",
    "where $\\vec u$ is the flux, $\\mathcal K$ the permeability tensor and $p$ the fulid pressure. Further, the conservation of mass gives\n",
    "$$\n",
    "\\frac{\\partial \\phi \\rho}{\\partial t} + \\nabla \\cdot \\rho \\vec u = q,\\quad \\text{in}\\ \\Omega \\\\\n",
    "u\\cdot n = 0,\\quad \\text{on}\\ \\partial \\Omega\n",
    "$$\n",
    "for porosity $\\phi$, fluid density $\\rho$, and source/sink term $q$.\n",
    "\n",
    "To solve this system of equation we need a constitutive law relating the fluid density to the pressure:\n",
    "$$\n",
    "\\rho = \\rho_r e^{c(p - p_r)},\n",
    "$$\n",
    "for reference density $\\rho_r$ and pressure $p_r$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import statements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.sparse as sps\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Porepy modules\n",
    "import porepy as pp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define constitutive laws and constants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We set the porosity to 0.2 and let set the permeability to the default value (i.e. $\\mathcal K = 1$).\n",
    "We define the depenecy of $\\rho$ on $p$ as a function. Note that we have to use the exponent function ad.exp (and not np.exp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define data\n",
    "dt = 0.2                           # Time step\n",
    "phi = 0.2                          # Porosity \n",
    "c = 1e-1                           # Compressibility\n",
    "\n",
    "# Constitutive law\n",
    "def rho(p):\n",
    "    rho0 = 1\n",
    "    p_ref = 1\n",
    "    return rho0 * pp.ad.functions.exp(c * (p - p_ref))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discretization\n",
    "\n",
    "We use a finite-volume method to discretize the model equation. As a first step we create a partition of the domain into grid cells:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVKklEQVR4nO3de5RlZXnn8e+PbrBpaAWCYXFTMGIcQsaYoKBkJQQcRSXCZGWIZlTMJOkkExRzWUaNiawYL5NJHF0rV4IoBgODDSiTyQVFjbkYoqATkTai3LqhgVYHwVsEfPLH3tU5FlXV1ZWu81bX+/2sdVafs/c++3neU9W/s897Tp2dqkKS1I+9WjcgSZoug1+SOmPwS1JnDH5J6ozBL0mdMfglqTMGv1aMJF9O8rh51r0kyd9OuydpNTL4O5Pkx5N8bAzZbUn+Isn3/zv2V0keP2vZhiRvTnJrkq8kuT3JpiQnLLSvqtq/qm5eYh/7JDkvyU1jzVuTXJjkqKXsb7mMPV7cuo/llOSo8fdibeteNDeDvyNJfhF4C/AG4BDgMcDvA2csYV9z/qdO8gjgA8B3A6cDjwT+A3Ap8Oxd2dcu2gQ8D/hx4FHAk4DrgFN3w761m/mk0FhVeengwhCGXwb+ywLbPBX4CHAvsA34XWCfifUF/DxwE3AL8OFx2VfGff8Y8FPjfffbST/fsq+JZY8fr38bcBVwH/CPwOuAv51nX88AvgYcuUC9w8b9fRH4LPDTE+vOA94NXAzcD3wSeALwKuAeYAvwzIntPwS8cezrPuC9wEHjupOBrbNq3zr2eBrwDeCB8fH6fxM/m7eNj9sdwG8Ca+YZx77ARcD/BzYDr5isN47zcmD7+DN62cS6RzA88d85Xt4CPGKy73F/94y9nAk8B/jM+Li9emJfewGvBD4HfAG4bOIxuH38WX55vDwNeAnwd8D/Grd/w7jP757Y57cDXwUe3fr/y2q/NG/Ay5R+0EPoPAisXWCb7wNOBNYCR43B8vKJ9QW8DzgI2Hdi2eMntrkUeMci+llwX+N+LgP2A44bA3G+4H8T8Nc7qfdhhlc364DvGYPxlHHdecDXgWeNY3/nGJq/CuwN/DTjk9O4/YfGfo4b+7scuHhcdzLzBP9ErYtnrb8S+KNxX9/O8ITyMwuNFTgQOAL4p5l6YxhfB/w6sA/wOOBm4Fnj+t8A/mGs8Wjg74HXTfT94HjfmTFvB/4U2AB8F8OT69Hj9ueO+zqC4Qnlj4BLxnVHjT/LtRN9v2Tc/0vHx3jf8efxPya2ORf4P63/r/Rwad6Alyn9oOG/Anft4n1eDlw5cbtmwnLWssngfz/wponb38PwCuI+4J8Xsy9gDcNR8RMn1r2B+YP/j4FLFxjHkcBDwIaJZW9kfIIaw/h9E+t+mOFIdc14e8PY2wHj7Q/NGuOxDEfya9jF4GeYcvsXxie/cdkLgA/OM5YdQT7e/in+LfhPAG6ftf2rgLeP1z8HPGdi3bOAW8frJzME++wxnzCx/XXAmeP1zcCpE+sOHX9mMwcNcwX/7N5OYHh1kPH2x4CzWv9f6eHiPFs/vgAcnGRtVT041wZJngC8GTgeWM/wn/i6WZttWUSdQ2duVNUngAOSPAO4YJH7evRYe3L9bTup+YQF1h8GfLGq7p+1v+Mnbt89cf1rwOer6qGJ2wD7MzyJMUdvewMHL9DDfB473ndbkpllezH/Y3PYrHWT1x8LHJbk3olla4C/mbjv5ON427hsxhfmGPPsx2X/iVpXJvnmxPqHGJ7I5vMtY6qqa5N8FTg5yTaGJ/2rFri/dhPf3O3HRxiOLM9cYJs/AD4NHFNVjwReDWTWNjv7OtdrgGcm2W8RPc23r+0M0wJHTix7zAL7eT/w1CRHzLP+TuCgJBtm7e+ORfQ4n9m9PQB8nuH9jvUzK5KsYXgimzF7zFsYfi4HV9UB4+WRVfVd89TdxjC9MlcfWximpA6YuGyoqueM6+9kCOzJvu9ccJTz2wI8e1atdVV1xxxjnDHX8ouAFwIvAjZV1deX2I92gcHfiar6EsP87e8lOTPJ+iR7J3l2kt8aN9vAMCXz5SRPBH5uEbu+m2EuecY7GcLpyiTHJVmTZB3fenS9s14fAq4Azhv7PBY4e4Ht38/wfsGVSb4vydrxI6U/m+S/VdUWhvnsNyZZl+Q/Aj/J8GbuUr0wybFJ1jPMnW8a+/4MsC7Jc5PsDbyGYQ58xt3AUUn2GnvfBlwN/E6SRybZK8l3JPnBeepeBrwqyYFJDgfOmVj3j8D9SX4lyb7jY39ckqeM6y8BXpPk0UkOZvh9WOpj8IfA65M8FmDc58ynw7YD3+Rbfy/mczHwnxnC/51L7EW7yODvSFX9DvCLDGG0neGo7RzgPeMmv8zwccj7GebN//cidnsecFGSe5OcNR6x/RBwI/B/Gef2gacAZ+1Cu+cwTCvcBbwDePtOtv9R4M/Hnr8E3MDwZPP+cf0LGOae72R4M/W14xPGUv3J2NddDG8Yvwx2PMH+d4ZprTsYXgFsnbjfu8d/v5Dk+vH6ixnejL2R4dM6m5iYLpvlN8b93TKObRPDK4aZJ8zTGd5XuYXhFcgFDJ8aguHTQh9jeEP4k8D147KleCvDtMzVSe5neKP3hLGPrwKvB/5u/L04cb6djE/K1zO8Gvib+bbT7jXzpoqkRUryIYY3aGe/Z9Gil58Dnl9V871CWPGSXAjcWVWvad1LL3xzV9qDJDmUYQrlI8AxwC8x/L3FHmn8y+ofAZ7cuJWuONUj7Vn2YfjM/P0MfyH9XobPw+9xkryOYUruf1bVLa376YlTPZLUGY/4Jakze8Qc/15Jk9clYecfWl9ttR3z6q/bsrZjnrqqqocd4O8RwV/09wPr8Ze0tzH7WPdRu/GYZ/8BJuBUjyR1x+CXpM4Y/JLUGYNfkjpj8EtSZwx+SeqMwS9JnTH4JakzBr8kdcbgl6TOGPyS1BmDX5I6Y/BLUmcMfknqjMEvSZ1ZtuBPcmGSe5LcMLHsoCTvS3LT+O+By1VfkjS35Tzifwdw2qxlrwSuqapjgGvG25KkKVq24K+qDwNfnLX4DOCi8fpFwJnLVV+SNLdpn3rxkKraNl6/Czhkvg2TbAQ27ri9zI3NZW2jui1rO+bVX7dlbce8MjQ7525VVZJ5T0VZVecD5wPEk62v+rota/dWt2Vtxzz92nOZ9qd67k5yKMD47z1Tri9J3Zt28F8FnD1ePxt475TrS1L3Uss0iZLkEuBk4GDgbuC1wHuAy4DHALcBZ1XV7DeA59qXUz2rvG7L2r3VbVnbMTeoXfWwGZ9lC/7dyeBf/XVb1u6tbsvajrlB7TmC37/claTOGPyS1BmDX5I6Y/BLUmcMfknqjMEvSZ0x+CWpMwa/JHXG4Jekzhj8ktQZg1+SOmPwS1JnDH5J6ozBL0mdaXbqxV0R+jtXZo/nB+1tzD7WfdT2nLtLVHT6Pdod1W1Zu7e6LWs75unXnotTPZLUGYNfkjpj8EtSZwx+SeqMwS9JnTH4JakzBr8kdcbgl6TOGPyS1BmDX5I6Y/BLUmcMfknqjMEvSZ0x+CWpMwa/JHWmSfAn+YUkn0pyQ5JLkqxr0Yck9WjqwZ/kcOBlwPFVdRywBnj+tPuQpF61mupZC+ybZC2wHrizUR+S1J2pn3qxqu5I8tvA7cDXgKur6urZ2yXZCGzccXt6Le7g+UH7qN1b3Za1HfPKkKrpng0yyYHA5cCPAfcC7wY2VdXFC9xnyl2OdfH8oD3U7q1uy9qOuUHtqoc977SY6nkGcEtVba+qB4ArgKc36EOSutQi+G8HTkyyPkmAU4HNDfqQpC5NPfir6lpgE3A98Mmxh/On3Yck9Wrqc/xL4Rz/6q/bsnZvdVvWdswNaq+QOX5JUkMGvyR1xuCXpM4Y/JLUGYNfkjpj8EtSZwx+SeqMwS9JnTH4JakzBr8kdcbgl6TOGPyS1BmDX5I6M/VTLy5F6O+UaT2eJq63MftY91F7JZ56cY8I/qLTr1PtqG7L2r3VbVnbMU+/9lyc6pGkzhj8ktQZg1+SOmPwS1JnDH5J6ozBL0mdMfglqTMGvyR1xuCXpM4Y/JLUGYNfkjpj8EtSZwx+SeqMwS9JnTH4JakzTYI/yQFJNiX5dJLNSZ7Wog9J6lGrE7G8FfjLqvrRJPsA6xv1IUndSdV0zw2T5FHAJ4DH1SKLJ5lyl2NdPFtQD7V7q9uytmNuULvqYSfianHEfzSwHXh7kicB1wHnVtVXJjdKshHYuOP2VFsceH7QPmr3Vrdlbce8MrQ44j8e+AfgpKq6Nslbgfuq6tcWuI9H/Ku8bsvavdVtWdsxN6g9xxF/izd3twJbq+ra8fYm4Hsb9CFJXZp68FfVXcCWJN85LjoVuHHafUhSr1p9quelwLvGT/TcDPxEoz4kqTtTn+NfCuf4V3/dlrV7q9uytmNuUHuFzPFLkhoy+CWpMwa/JHVmp8Gf5KVJDpxGM5Kk5beYI/5DgI8muSzJaUlW2h+hSZJ2waI+1TOG/TMZPnZ5PHAZ8Laq+tzytrejvp/qWeV1W9burW7L2o65Qe2lfqpn/DK1u8bLg8CBwKYkv7U7m5QkLb+dHvEnORd4MfB54ALgPVX1QJK9gJuq6juWvUmP+Fd93Za1e6vbsrZjblB7id/OeRDwI1V12+TCqvpmktN3U3+SpCnxL3cXqotHJz3U7q1uy9qOuUFt/3JXkmTwS1JnDH5J6ozBL0mdafV9/Lsk9HeuzB7PD9rbmH2s+6i9Es+5u0cEf9Hpu/Ed1W1Zu7e6LWs75unXnotTPZLUGYNfkjpj8EtSZwx+SeqMwS9JnTH4JakzBr8kdcbgl6TOGPyS1BmDX5I6Y/BLUmcMfknqjMEvSZ0x+CWpMwa/JHWmWfAnWZPk40n+rFUPktSjlkf85wKbG9aXpC41Cf4kRwDPBS5oUV+Setbq1ItvAV4BbJhvgyQbgY07bi9/Tw/j+UH7qN1b3Za1HfPKMPXgT3I6cE9VXZfk5Pm2q6rzgfPH+1SP58p0zNZdbbUd8/Rrz6XFVM9JwPOS3ApcCpyS5OIGfUhSl1JNjqXH4sMR/y9X1ek72c4j/lVet2Xt3uq2rO2YG9SuetiBv5/jl6TOND3iXyyP+Fd/3Za1e6vbsrZjblDbI35JksEvSZ0x+CWpMwa/JHXG4Jekzhj8ktQZg1+SOmPwS1JnDH5J6ozBL0mdMfglqTMGvyR1xuCXpM60OvXiLgn9nTKtx9PE9TZmH+s+anvqxSUqOv061Y7qtqzdW92WtR3z9GvPxakeSeqMwS9JnTH4JakzBr8kdcbgl6TOGPyS1BmDX5I6Y/BLUmcMfknqjMEvSZ0x+CWpMwa/JHXG4Jekzhj8ktQZg1+SOjP14E9yZJIPJrkxyaeSnDvtHiSpZ6ma7ikCkhwKHFpV1yfZAFwHnFlVNy5wnyl3OdbFk0b0ULu3ui1rO+YGtasedj6WqR/xV9W2qrp+vH4/sBk4fNp9SFKvmp56MclRwJOBa+dYtxHYuOP29NrawfOD9lG7t7otazvmlWHqUz07Cif7A38NvL6qrtjJtk71rPK6LWv3VrdlbcfcoPZKmOoBSLI3cDnwrp2FviRp92rxqZ4AbwM2V9Wbp11fknrX4oj/JOBFwClJPjFentOgD0nqUrM5/l3hHP/qr9uydm91W9Z2zA1qr5Q5fklSOwa/JHXG4Jekzhj8ktQZg1+SOmPwS1JnDH5J6ozBL0mdMfglqTMGvyR1xuCXpM4Y/JLUGYNfkjpj8EtSZ5qec3exQn/nyuzx/KC9jdnHuo/aK/Gcu3tE8Bedfo92R3Vb1u6tbsvajnn6tefiVI8kdcbgl6TOGPyS1BmDX5I6Y/BLUmcMfknqjMEvSZ0x+CWpMwa/JHXG4Jekzhj8ktQZg1+SOmPwS1JnDH5J6ozBL0mdaRL8SU5L8s9JPpvklS16kKRepWq6pwhIsgb4DPCfgK3AR4EXVNWNC9xnyl2OdfGkET3U7q1uy9qOuUHtqoedj6XFEf9Tgc9W1c1V9Q3gUuCMBn1IUpdanHrxcGDLxO2twAmzN0qyEdi44/by9zWnlufKdMzWXY21HfNUzfliY8Wec7eqzgfOB0jysao6vnFLU+WYV7/exguOeaVoMdVzB3DkxO0jxmWSpCloEfwfBY5JcnSSfYDnA1c16EOSujT1qZ6qejDJOcBfAWuAC6vqUzu52/nL39mK45hXv97GC455RZj6xzklSW35l7uS1BmDX5I6s6KDv7evdkhyZJIPJrkxyaeSnNu6p2lJsibJx5P8WetepiHJAUk2Jfl0ks1Jnta6p+WW5BfG3+sbklySZF3rnna3JBcmuSfJDRPLDkryviQ3jf8e2LJHWMHBP361w+8BzwaOBV6Q5Ni2XS27B4FfqqpjgROBn+9gzDPOBTa3bmKK3gr8ZVU9EXgSq3zsSQ4HXgYcX1XHMXyw4/ltu1oW7wBOm7XslcA1VXUMcM14u6kVG/x0+NUOVbWtqq4fr9/PEAaHt+1q+SU5AngucEHrXqYhyaOAHwDeBlBV36iqe5s2NR1rgX2TrAXWA3c27me3q6oPA1+ctfgM4KLx+kXAmdPsaS4rOfjn+mqHVR+CM5IcBTwZuLZxK9PwFuAVwDcb9zEtRwPbgbeP01sXJNmvdVPLqaruAH4buB3YBnypqq5u29XUHFJV28brdwGHtGwGVnbwdyvJ/sDlwMur6r7W/SynJKcD91TVda17maK1wPcCf1BVTwa+wgp4+b+cxnntMxie9A4D9kvywrZdTV8Nn59v/hn6lRz8XX61Q5K9GUL/XVV1Ret+puAk4HlJbmWYzjslycVtW1p2W4GtVTXzam4TwxPBavYM4Jaq2l5VDwBXAE9v3NO03J3kUIDx33sa97Oig7+7r3ZIEoZ5381V9ebW/UxDVb2qqo6oqqMYfsYfqKpVfSRYVXcBW5J857joVGDe81GsErcDJyZZP/6en8oqf0N7wlXA2eP1s4H3NuwFWNnfzrmUr3bY050EvAj4ZJJPjMteXVV/3q4lLZOXAu8aD2puBn6icT/LqqquTbIJuJ7h02sfZwV+lcG/V5JLgJOBg5NsBV4LvAm4LMlPArcBZ7XrcOBXNkhSZ1byVI8kaRkY/JLUGYNfkjpj8EtSZwx+SeqMwS9JnTH4JakzBr+0BEmekuSfkqxLst/4PfPHte5LWgz/gEtaoiS/CawD9mX47p03Nm5JWhSDX1qi8esWPgp8HXh6VT3UuCVpUZzqkZbu24D9gQ0MR/7SHsEjfmmJklzF8FXSRwOHVtU5jVuSFmXFfjuntJIleTHwQFX96Xh+6L9PckpVfaB1b9LOeMQvSZ1xjl+SOmPwS1JnDH5J6ozBL0mdMfglqTMGvyR1xuCXpM78K1UoeLMDmcRNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Create grid\n",
    "g = pp.CartGrid([11,11])\n",
    "g.compute_geometry()\n",
    "pp.plot_grid(g, plot_2d=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, the model equation is integrated over each controll volume (i.e., each cell of the grid), and the divergence theorem is applied to the flux term:\n",
    "$$\n",
    "\\int_\\Omega \\phi \\frac{\\partial \\rho}{\\partial t} dV - \\int_{\\partial\\Omega}\\vec n\\cdot(\\rho\\vec u)dS - \\int_\\Omega q dV= 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The key-point of the finite-volume discretization is how the flux-term $\\vec u$ is approximated. We do not cover that in this tutorial(see e.g., I. Aavatsmark. An introduction to multipoint flux approximations for quadrilateral grids. Comput. Geosci., Vol. 6, No. 3, pp. 405–432, 2002. DOI: 10.1023/A:1021291114475).\n",
    "However, the main idea is that the fluid flux $\\vec u$ across a face is expressed as a linear combination of the cell-centered pressures $\\vec u = \\text{flux}\\ \\vec p$. Here, $\\text{flux}$ is the discretization matrix and $\\vec p$ is the vector of all cell-centered pressures.\n",
    "\n",
    "In porepy we can obtain the discretization matrix with, e.g, the two-point flux approximation:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize default data (i.e., unit parameters)\n",
    "data = pp.initialize_default_data(g, {}, 'flow')\n",
    "# Define flux discretization:\n",
    "flx_disc = pp.Tpfa('flow')\n",
    "# Discretize\n",
    "flx_disc.discretize(g, data)\n",
    "# The flux discretization can now be found in the dictionary as:\n",
    "flux = data[pp.DISCRETIZATION_MATRICES]['flow']['flux']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the negative sign in front of the surface-integral is included into the flux discretization matrix.\n",
    "\n",
    "The density is defined at the cell centers, but in the flux term we need to evaluate it at the faces. To do so, we simply take the average of the two neighbooring cells (note that other alternatives, such as upstream weighting, are commonly used).\n",
    "\n",
    "We also create discretized versions of the divergence operator div. The discrete divergence operator sums the fluxes in and out of each grid cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "cell_faces_T = g.cell_faces.T\n",
    "def div(x):\n",
    "    \"\"\"\n",
    "    Discrete divergence\n",
    "    \"\"\"\n",
    "    return cell_faces_T * x\n",
    "\n",
    "def avg(x):\n",
    "    \"\"\"\n",
    "    Averageing. Note that this is not strictly correct for the boundary faces since\n",
    "    these only have 1 cell neighboor, but we have zero flux condition on these, so \n",
    "    this is not a problem.\n",
    "    \"\"\"\n",
    "    return 0.5 * np.abs(g.cell_faces) * x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Residual function\n",
    "To discretize the time-deriveative, we use backward Euler. Further, we assume that the densities are constant over each cell so we can take them out of the integral:\n",
    "$$\n",
    "\\int_\\Omega \\phi \\frac{\\rho^k - \\rho^{k-1}}{\\Delta t} dV =\\phi \\frac{\\rho^k - \\rho^{k-1}}{\\Delta t} \\int_\\Omega dV = \\phi \\frac{\\rho^k - \\rho^{k-1}}{\\Delta t}V,\n",
    "$$\n",
    "where $V$ is the volume of the cell. The same is also done for the source term.\n",
    "\n",
    "This gives us the residual\n",
    "$$\n",
    "\\phi \\frac{\\rho^k - \\rho^{k-1}}{\\Delta t} V + \\text{div}(\\text{avg}(\\rho^k)\\text{flux } p^k) - q^k V= 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(p, p0):\n",
    "    # darcy:\n",
    "    u = flux * p\n",
    "\n",
    "    # Source:\n",
    "    src = np.zeros(g.num_cells)\n",
    "    src[60] = 1\n",
    "\n",
    "    # Define residual function\n",
    "    time = phi * (rho(p) - rho(p0)) / dt * g.cell_volumes\n",
    "    advection = div(avg(rho(p)) * u)\n",
    "    lhs = time + advection\n",
    "    rhs = src * g.cell_volumes\n",
    "\n",
    "    return lhs - rhs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initialize Ad variable\n",
    "To initialize an AD variable create an Ad_array(...) with values equal the initial value and jacobian equal the identity matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set initial condition\n",
    "p0 = np.zeros(g.num_cells)\n",
    "p = pp.ad.Ad_array(p0, sps.diags(np.ones(p0.shape)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time loop\n",
    "We are now ready to set up the time loop. We will set up a simple Newton iteration to find the zero of the residual function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solving time step:  1\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/eke001/Dropbox/workspace/python/ppdir/src/porepy/viz/plot_grid.py:241: UserWarning: Attempting to set identical bottom == top == 0.0 results in singular transformations; automatically expanding.\n",
      "  ax.set_zlim3d(z)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARoAAAD9CAYAAABjhbg0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAo+ElEQVR4nO3df3gU1b348fcnu1kCAaQR0JiAgET5ESw/guitT78CopS2sQq1ULAq5HqvX73Fr5f20mot2F6lWitYvNdSFVEficptG1TECtYrohEjqVXTImgCSUCQnwJJSHZzvn/MLC6BZHeSnezO8nk9zz5kd85+5swO+9kzZ+bMEWMMSinlprREV0Aplfo00SilXKeJRinlOk00SinXaaJRSrlOE41SynWaaNogIkdEZFAry24QkTc7u05KeVHCE42IfF9Eyuwv9S4ReVlELu1APCMig1u81kNEfiMiVSJyVER2iMgqERnXVixjTHdjzKftrEdARBaIyFZ7nVUi8riIDGhPPLfYdXw60fVwk4gMsP9f+BNdl9NVQhONiNwOLAbuAc4C+gP/BVzVjlin/E8kIl2A14ARwLeAnsBQoBj4hpNYDq0CCoHvA2cAXwXeAybGIbaKM01CLjPGJOSB9eU7Any3jTIXAW8DB4FdwFIgELHcALcAW4FK4A37taN27O8BRfZ7M6PU54RYEa8Ntv8+E1gNfAFsAn4BvNlKrMuBeqBfG+s7x463H9gG/HPEsgXA88DTwGHgA+B84CfAHqAauCKi/OvAvXa9vgBKgCx72WVATYt1V9l1nAw0Ak325/V+xL55zP7caoFfAr5WtqMrsAI4APwd+HHk+uzt/B/gc3sf/TBiWResH5qd9mMx0CWy3na8PXZdvgNMAT62P7efRsRKA+YDnwD7gOciPoMd9r48Yj8uAW4ANgIP2uXvsWOOiIjZF6gD+iTqe5Iqj8St2PpPHgT8bZQZA1wM+IEB9n/k2yKWG+BVIAvoGvHa4IgyxcATMdSnzVh2nOeATCDf/gK2lmgWAf8bZX1vYLXeMoCR9hdxgr1sAdAAXGlv+5P2l/QOIB34Z+xkaJd/3a5Pvl2//wGetpddRiuJJmJdT7dY/kfgd3asvlgJ7F/a2lbgK0Au8Lfw+uwv/3vAXUAAGAR8ClxpL78bKLXX0Qd4C/hFRL2D9nvD2/w58AzQAxiOlcwH2uXn2rFysRLY74CV9rIB9r70R9T7Bjv+v9mfcVd7f/wqosxc4IVEf0lT4ZG4FcNM4DOH77kN+GPEcxP+crZ4LTLRrAMWRTwfidVC+gLYEksswIf1qz8kYtk9tJ5ofg8Ut7Ed/YAQ0CPitXuxE6L95X81Ytm3sX6JffbzHnbdetnPX2+xjcOwWio+HCYarEPYY9jJ1n5tBvCXVrbleOKwnxfxZaIZB+xoUf4nwHL770+AKRHLrgSq7L8vw0okLbd5XET594Dv2H//HZgYsSzb3mfhH6lTJZqWdRuH1foR+3kZcG2iviOp9Ejkcek+oLeI+I0xwVMVEJHzgd8ABUA3rP8077UoVh3DerLDT4wxfwV6icjlwKMxxupjrzty+fYo6zy/jeXnAPuNMYdbxCuIeL474u96YK8xJhTxHKA7VtLkFHVLB3q3UYfWnGu/d5eIhF9Lo/XP5pwWyyL/Phc4R0QORrzmAzZEvDfyc9xuvxa27xTb3PJz6R6xrj+KSHPE8hBW4mzNCdtkjHlHROqAy0RkF9aPzOo23q9ilMjO4Lexfjm/00aZ/wb+AeQZY3oCPwWkRZlow8/XA1eISGYMdWot1udYzex+Ea/1byPOOuAiEcltZflOIEtEerSIVxtDHVvTsm5NwF6s/qpu4QUi4sNKnGEtt7kaa7/0Nsb0sh89jTHDW1nvLqzDlVPVoxrrEK9XxKOHMWaKvXwnVoKIrPfONreyddXAN1qsK8MYU3uKbQw71esrgFnAdcAqY0xDO+ujIiQs0RhjDmEdfz8sIt8RkW4iki4i3xCR++xiPbAOcY6IyBDg5hhC78bqCwh7EuvL8EcRyRcRn4hkcGLrIVpdQ8AfgAV2PYcB17dRfh1Wf88fRWSMiPjtU+z/KiKzjTHVWP0R94pIhohcCMzB6vxtr1kiMkxEumH1fayy6/0xkCEi3xSRdOBOrD6MsN3AABFJs+u+C/gz8ICI9BSRNBE5T0T+TyvrfQ74iYh8RURygFsjlm0CDovIf4hIV/uzzxeRsfbylcCdItJHRHpj/X9o72fwCPCfInIugB0zfPbyc6CZE/9ftOZp4GqsZPNkO+uiWkjo6W1jzAPA7Vj/+T/H+lW6FfiTXWQe1unhw1j9Hs/GEHYBsEJEDorItfYv0nigAngJu28GGAtc66C6t2I10z8DngCWRyk/DVhj1/kQ8CFWcltnL5+B1XewE6vz9ed2gmqvp+x6fYbVwfxDOJ7Q/y/WYWItVgunJuJ9z9v/7hORzfbfP8DqvK3AOpu0iojDzxbutuNV2tu2CqtFFE7Q38LqF6vEamE9inVWC6yzWWVYHcgfAJvt19pjCdZhzp9F5DBWx/A4ux51wH8CG+3/Fxe3FsT+EdiM1drZ0Fo55Uy400t5mIi8jtWh27LPKRF1uRmYboxprQWU9ETkcWCnMebORNclVehFSqpDRCQb65DkbSAP+Hes6508yb5y+xpgVIKrklISPgRBeV4A65qVw1hXYJdgXY/iOSLyC6xD3PuNMZWJrk8i2MNk9ojIh60sFxF5SES2icjfRGR0THH10EkpFSYiX8e6ZutJY0z+KZZPwbrIcQpWH9gSY0ybYwZBWzRKqQjGmDewhmK05iqsJGSMMaVY16S1dqLguGh9NNrcUcp9La8Nc2SwiKmLsewu+AhreEvYMmPMMgery+HECx1r7Nd2tfUm7QxWyuPqgH+JsewCaDDGxHwNWbxoolHK44RO/SLXcuLV37nEcEW79tEo5XFpWEPPY3nEwWrgB/bZp4uBQ/bV5G3SFo1SHidYo2DjEktkJdbI+d4iUgP8PBzeGPMI1tXuU7DuoVQH3BhLXE00SnlcPA+djDEzoiwP3yDOEU00SnlcPFs0btFEo5THdXJncLske/2UUlF4oUWjZ51SzLvvvsuFF15IQ0MDR48eZfjw4Xz44SmHragU0clnndpFWzQpZuzYsRQWFnLnnXdSX1/PrFmzyM8/aciKSiFeaNFEG1SpQxA8qLGxkbFjx5KRkcFbb72Fz+dLdJVU2zo0BCFPxDwUY9kp8J5eGaziYt++fRw5coSmpiYaGhrIzIzldsnKq7RFoxKisLCQ6dOnU1lZya5du1i61LP3oTpddKhFc4GI+V2MZcdri0bFw5NPPkl6ejrf//73CYVC/NM//ROvvfYaEyZMSHTVlEvCncHJTFs0SiVeh1o0w0RMrFNHjNEWjVKqPfSCPaWU67zQGayJRimP0xaNUsp12qJRSrlOSP6zTppolPI4AdJj/SYH3axJ6zTRKOVxIuDXRKOUcpMIpCf5cDZNNEp5nKMWTYIkefWUUtGIQHqXRNeibZpolPI6D1xIk7J32Fu7di0XXHABgwcPZtGiRXGLW11dzfjx4xk2bBjDhw9nyZIlcYsdFgqFGDVqFN/61rfiGvfgwYNMmzaNIUOGMHToUN5+++24xX7wwQcZPnw4+fn5zJgxg4aGhuhvasPs2bPp27fvCTft2r9/P5MmTSIvL49JkyZx4MCBjlY7NYQTTSyPBEnJRBMKhbjlllt4+eWXqaioYOXKlVRUVMQltt/v54EHHqCiooLS0lIefvjhuMUOW7JkCUOHDo1rTIC5c+cyefJk/vGPf/D+++/HbR21tbU89NBDlJWV8eGHHxIKhSguLu5QzBtuuIG1a9ee8NqiRYuYOHEiW7duZeLEiXH9AfE8TTSdb9OmTQwePJhBgwYRCASYPn06JSUlcYmdnZ3N6NGjAejRowdDhw6ltjbqjKAxq6mp4aWXXqKoqChuMQEOHTrEG2+8wZw5cwAIBAL06tUrbvGDwSD19fUEg0Hq6uo455xzOhTv61//OllZWSe8VlJSwvXXXw/A9ddfz5/+9KcOrSNlCOCL8ZEgKZloamtr6dfvy+mBc3Nz45oMwqqqqigvL2fcuHFxi3nbbbdx3333kZYW311TWVlJnz59uPHGGxk1ahRFRUUcPXo0LrFzcnKYN28e/fv3Jzs7mzPOOIMrrrgiLrEj7d69m+zsbADOPvtsdu/eHfd1eJIeOqWuI0eOMHXqVBYvXkzPnj3jEvPFF1+kb9++jBkzJi7xIgWDQTZv3szNN99MeXk5mZmZcTv0OHDgACUlJVRWVrJz506OHj3K00/HeoeU9hERRDp0G5fUIUCXGB8J4slE09zcTFs37MrJyaG6uvr485qaGnJycuK2/qamJqZOncrMmTO55ppr4hZ348aNrF69mgEDBjB9+nRee+01Zs2aFZfYubm55ObmHm99TZs2jc2bN8cl9rp16xg4cCB9+vQhPT2da665hrfeeisusSOdddZZ7NplzSe/a9cu+vbtG/d1eJK2aNzR2NhIY2Njq8lm7NixbN26lcrKShobGykuLqawsDAu6zbGMGfOHIYOHcrtt98el5hh9957LzU1NVRVVVFcXMyECRPi1jI4++yz6devH1u2bAFg/fr1DBs2LC6x+/fvT2lpKXV1dRhjWL9+vSud2YWFhaxYsQKAFStWcNVVV8V9HZ7kgUST5GffW7djxw769+9PIBA4qQnt9/tZunQpV155JaFQiNmzZzN8+PC4rHfjxo089dRTjBgxgpEjRwJwzz33MGXKlLjEd9Nvf/tbZs6cSWNjI4MGDWL58uVxiTtu3DimTZvG6NGj8fv9jBo1iptuuqlDMWfMmMHrr7/O3r17yc3NZeHChcyfP59rr72Wxx57jHPPPZfnnnsuLvVPCUk+BMGT9wxuaGjg7bff5uKLL6a2tpbzzjtPj9eVl3XoP29BDzFlMXbryf/qPYMdS0tLO94Xo8lGnbb0yuDOUV1dzSeffHJSn83kyZNdW6dXY7sd36uxPc0DZ52SPA/G7lQtm71797q2Pq/Gdju+V2N7mgdaNElePWeqq6sZMmQ4oVDj8deiH075gJCDtXxZPrZDNafxrfc4Owx0vg1uxhfxO6xP7Ov4st5Ot9mPMU0O6+QRmmg6n5Vk1jh4xxTgVQflJwF/cVB+PPCmg/IAlwKlDspfDLznoPyYdsR3sg2X4uwzAutzcrofnO7nFBUegpDEUqKPRqnTWpyvoxGRySKyRUS2icj8UyzvLyJ/EZFyEfmbiETN4inXolHqtBPuDI5HKBEf8DBWk7EGeFdEVhtjIm9RcCfwnDHmv0VkGFbTckBbcbVFo5TXxbdFcxGwzRjzqTGmESgGWl6CbYDwAL8zgJ3RgmqLRimvc9YZ3FtEyiKeLzPGLIt4ngNURzyvAVrenmAB8GcR+TcgE7g82ko10SiVCmL/Ju+Nw5XBM4AnjDEPiMglwFMikm+Mae549ZRSySm+Z51qgX4Rz3Pt1yLNASYDGGPeFpEMoDewp7Wg2kejlNfFt4/mXSBPRAaKSACYDqxuUWYHMBFARIYCGcDnbQX1ZItm3759BIMJmnJPqWQTx7NOxpigiNwKvILVTnrcGPORiNwNlBljVgP/DvxeRP4fVsfwDSbK6GxPJppu3bpx7Ngxtm/fnuiqKJV4cb4y2BizhhZXQxpj7or4uwL4mpOYnkw0Xbt2pVu3bhw+fJj6+noyMjJ05LY6fekQBPeICPn5+ezdu5e6ujq6du1q39Dbh7PLzX1Y1yY5KT/eYflLHZQPv+dih+Wd3Ge4PfGdbIPTzyj8Hqf7wel+TlGaaNwXCATw+XzU1dXRpUsXrIF2yTR26VKcjUMCK2k4mStqGPCxg/LntyO+07FU7Rnf5XQ/ON3PKSzJ86jnEw2Az+cjMzOT+vr6RFdFqc6nLZrOIyJ069Yt0dVQSaygoIDevXufNAOm58XxrJNb9DoaddooKytrNcl4eq5vD8yCoIlGKTw+17cmGqW8wdNzfXtg7u2U6aNRKt48M9e3dgYrlRqSeq5vwRptlMT00EmpVnhmrm8PHDppolGqFZ6Z61s7g5XyhhkzZnDJJZewZcsWcnNzeeyxx5g/fz6vvvoqeXl5rFu3jvnzT7pPd/JI8kSTgn00fpJr7JLTcUjh9wxzWP58l+M7HUvVnvFdTveD0/3cupUrV57y9fXr1ztYR4J4YLqVFEw0QZJr7NIYoNJBeYCBwD4H5c8E6hyU79aO+E62YSDtG9/l5txRTgd5eoiedVJKuc4DQxA00SjlddqiUUq5ThONUsp1mmiUUp1CzzoppVylLRqllOv0rJNSynXaonFHKBRKdBWUSh6aaNyxf/9+jhw5QlVVFcaY5B2+r1Rn0ETjjj59+pCZmQnA0aNH8fv9BAIBe14npU4/Rs86uUNEGDBgADt37qSpqYn6+nq7ZeMnuQZJ+rHG/jjhxxpf5KS8kxkg2hPfyTa0dyCpm5PUJfk3sQNMGjQm+Y2vPJtoIqWnp5Oenm733QRJrkGSA3E24BGgG/jbnDP9REGBbAfld4nz+I4HbbZnIKmbk9RdmrLTrRiBoC/W1nyzq3VpTUokmjCfL3V/tVTHlZWVxVTuwQcf5NFHH0VEGDFiBMuXLycjI3mbDEaEkD/Wr3Kjq3VpjXZqKBWhtraWhx56iLKyMj788ENCoRDFxcWJrlZUIZ8vpkeipFSLRql4CAaD1NfXk56eTl1dHeecc06iq9QmgxBK8j4obdEoFSEnJ4d58+bRv39/srOzOeOMM7jiiisSXa02GYQgvpgeiaKJRqkIBw4coKSkhMrKSnbu3MnRo0d5+umnE12tNhmERrrE9EgUTTRKRVi3bh0DBw6kT58+pKenc8011/DWW28lulptCh86xfJIFE00SkXo378/paWl1NXVYYxh/fr1DB06NNHViiqeiUZEJovIFhHZJiKnnPpBRK4VkQoR+UhEnokWUzuDlYowbtw4pk2bxujRo/H7/YwaNYqbbrop0dVqU7iPJh5ExAc8jDXFRA3wroisNsZURJTJA34CfM0Yc0BEos6sp4lGqRYWLlzIwoULE12NmFmHTnH7Kl8EbDPGfAogIsXAVUBFRJl/Bh42xhwAMMbsiRZUE41SHmd1BgdiLd5bRCKvXFxmjFkW8TwHqI54XgOMaxHjfAAR2Yg1tmOBMabNy61TMNH4Sa6xS07HIdnvCToZke63hhW4Gd/xWCqn47vcnqQuua8z6QgDTg6d9hpjCjq4Sj+QB1wG5AJviMgIY8zBtt6QYoJAqYPyF3NiqzCaYTiefM3JuCKwksB4B+/5i8BCB+V/3o74jsdGOfmMwBrk6XQ/ON3PqSquh061QL+I57n2a5FqgHeMMU1ApYh8jJV43m0tqJ51Usrj4nx6+10gT0QGikgAmA6sblHmT1itGUSkN9ah1KdtBU3BFo1Sp594XSNjjAmKyK3AK1jHm48bYz4SkbuBMmPManvZFSJSAYSAHxlj2mzCaqJRyuPiPdbJGLMGWNPitbsi/jbA7fYjJppolPI4g3AsyadB0ESjlMd5YfS2JhqlPE4TjVKqUyTyFhCx0ESjlMfFeQiCK5K7dkqpqPTQySUHDhygvr6ebdu20djYiIiQlpamE8mp05J11inmsU4J4clEk5mZSXp6Oj169ACgubmZUChEc3Mz1iY5udzch3U5e6zaMSeSo3FFgPity/5jlea3hhW4Fb9dY6OcfEbgfD/4cLqfU3a6FT10ckcgEMDv93PWWWfxySeftFgaxPn8QB87KH8+juc4cjLnElgDJB2OXVpg/iPm4gvkV87HRjmdN6o9c1k53g/O9nOs060cPHiQoqIiPvzwQ0SExx9/nEsuucTBujqfHjop5TFz585l8uTJrFq1isbGRurqnCbNzqV9NEp5zKFDh3jjjTd44oknAKv1HAgkd/+HFxKNjt5WKkJlZSV9+vThxhtvZNSoURQVFXH06NGY3nvXXXexePHi48/vuOMOlixZ4lJNvxQeghDLI1E00SgVIRgMsnnzZm6++WbKy8vJzMxk0aJFMb139uzZPPnkk4B1gqK4uJhZs2a5WV1AZ0FQynNyc3PJzc1l3Djr7pXTpk1j8+bNMb13wIABnHnmmZSXl/PnP/+ZUaNGceaZTs++tU+yJxrto1Eqwtlnn02/fv3YsmULF1xwAevXr2fYsNhPuxcVFfHEE0/w2WefMXv2bBdr+qV4zoLgFk00SrXw29/+lpkzZ9LY2MigQYNYvnx5zO+9+uqrueuuu2hqauKZZ6JOdxQXeh2NUh40cuTImK+5aSkQCDB+/Hh69eqFz9d5rYxkP+ukiUapOGpubqa0tJTnn3++09bpcLqVhNDOYKXipKKigsGDBzNx4kTy8vI6bb3hPppYHomSgi0aP87nBzrfYXyHcxw5mnMJx2OX0vxp1rACl+K3a94ox3NZOd0P7ZkHyl3Dhg3j00/bnAzAFdpHkxCn37xOzTqvUwxSeV4n7aNRSrnMC0MQNNEo5XF6HY1SynXWWSedbkUp5SI9dFJKdQpNNEopV2kfjVLKdXodjVLKdV4YgqCJRimP00OnTtLc3EwwGCQUCiW6KkolhB46uaCpqYmmpiY++ugjjhw5QlpaGj6fz76JtM7rFPf4p+G8TqFQiIKCAnJycnjxxRcdrKfz6eltlxw+fJjm5mays7M5dOhQi6VB4E0H0S7F+TxQlQ7KD8TxHEemm/OxRU7nXXI8dsnhXFaOPiOwPien+8HZfnZyj5klS5YwdOhQvvjiCwfrSAwvJBpP3iYiKyuLLl26kJWVleiqqBRUU1PDSy+9RFFRUaKrEjO9TYRSHnPbbbdx3333cfjw4URXJSbNpCX9EARPtmiUcsuLL75I3759GTPGyb1uEi+esyCIyGQR2SIi20RkfhvlpoqIEZGCaDG1RaNUhI0bN7J69WrWrFlDQ0MDX3zxBbNmzeLpp59OdNVaFc8+GhHxAQ8Dk4Aa4F0RWW2MqWhRrgcwF3gnlrjaolEqwr333ktNTQ1VVVUUFxczYcKEpE4yAIa49tFcBGwzxnxqjGkEioGrTlHuF8CvgIZYgmqiUcrzrCEIsTyA3iJSFvG4qUWwHKA64nmN/dqXaxMZDfQzxrwUaw310EmpVlx22WVcdtllia5GVA4PnfYaY6L2qbRGRNKA3wA3OHmfJhqlPM4gHIvfWKdaoF/E81z7tbAeQD7wuogAnA2sFpFCY0yrFyppolHK4+I8evtdIE9EBmIlmOnA94+vy5hDQO/wcxF5HZjXVpIBTTRKpYR4nXUyxgRF5FbgFaxxHo8bYz4SkbuBMmPM6vbE1USjlMfFewiCMWYNsKbFa3e1UvayWGKmYKLxY41fipXTicj8WONynJR3OplaOwYxOp3gzfEgSYeT5jn6jKB9E8I53c+pySCEmpN7+1Iw0QSBvzgoP57kGoQJ1pfU4SR1jgc9Oo3vdCCpk88I2jNI0vl+Tk2mWTjWkNxDEFIw0Sh1ejFGCAW1RaOUcpNBE41Syl3GCMEmTTRKKVcJzaHk/iond+2UUtEZQA+dlFKuahZoSO6vcnLXTikVm2CiK9C2lEo0Ot2KOi1ZN6RJailxPxpjDA0NDTQ0xHQPHnWaKigoYPLkyYmuRvyFE00sjwTxdIvGGENTUxPHjh2jS5cudOmS3FdHqsSKZbqV6upqfvCDH7B7925EhJtuuom5c+d2Qu06wABNia5E2zybaJqbmykrKyMYDJKZmYl9bwysTXJyuXl7xsy4OTYq/B6Hk9Q5HovkNL6TbXD6GYXf43Q/ON3P0fn9fh544AFGjx7N4cOHGTNmDJMmTWLYMCeT23UyAxxLdCXa5slEs2fPHurr68nPz6eioqLF0iDwqoNok0iusVFgfUlbbldbhgEfOyh/fjviuzm5G7Rv7JLT/RxddnY22dnZAPTo0YOhQ4dSW1ub/IkmyftoPJlosrKyyMzMpFevXomuikphVVVVlJeXM27cuERXpW2aaNzh93uy2spDjhw5wtSpU1m8eDE9e/ZMdHXapolGKe9pampi6tSpzJw5k2uuuSbR1YlOE41S3mKMYc6cOQwdOpTbb7890dWJXZInmpS4jkapeNm4cSNPPfUUr732GiNHjmTkyJGsWbMm+hsTqRlrGrdYHgmiLRqlIlx66aUYYxJdDWf00Ekp5TpNNEop12miUUp1Ck00SilXaYsmEfzEerm5pT1jZtwcGxV+j5NL3n1YwwrcjO/mnEvh9zjdD073c4pqBuoTXYm2pWCiCdJikr0oppBcY6PA+pKWOih/Mc7HIjmN7+acS9C+sUtO93OKMkCS34opBRONUqchPXRSSrlK+2iUUq7TRKOUcl14CEIS00SjVCrQFo1SylV66KSUcp0Hbk6eMreJMMZQV1eX6GqoJJbS062EYnwkSEokmnCS0Vt8qraUlZWxdu3aqOXWrl3LBRdcwODBg1m0aFEn1KyD4jyvk4hMFpEtIrJNROafYvntIlIhIn8TkfUicm60mJ5PNOEkEwgESE9PT3R1lMeFQiFuueUWXn75ZSoqKli5cuUpZtpIMgZrCEIsjyhExAc8DHwDa5zKDBFpOV6lHCgwxlwIrALuixbX002AxsZG6urq6NKlC36/H2MM3bufwZEjTi43dzo2yum8UX6cj/vxY13276S807mmnMZ3sg1OP6Pwe5yOXXK2nw8fPkyPHj3aLLVp0yYGDx7MoEGDAJg+fTolJSXJP91K/A6LLgK2GWM+BRCRYuAqIubnMcZEji8pBWZFC+rZRBOeQC4yyfh8Platepampia6du0aMalc/DQ1NREKhcjIyIh7bLCSp4i41jpraGggPT0dn8+dQYZ1dXVkZGSQlhb/xnK49dq1a1dH8efPn8++ffsYPHgwOTk59O3bt9VDqNraWvr163f8eW5uLu+8806H6+4qZ2edeotI5JSdy4wxyyKe5wDVEc9rgLbmm5kDvBxtpZ48dKqvr6euro4hQ4ackGSOHTvmapIxxtDY2Ojq1LvBYNC1JAAgIq7eqjIQCNDY2OhKbBEhIyOD+vp6R9uwaNEifv/735OdnU1tbS179uxJrU5hZ300e40xBRGPZaeMGQMRmQUUAPdHK+vJRHP48GG6du1KVlbWCUnm2LFjriUZgGPHjhEIBFyLb4zBGONKayBMRGhubnYtvt/vJxQKuZbMfD4fgUCAhgbnl8L+5je/iZpscnJyqK7+8ge9pqaGnJycDtXZdeHT27E8oqsF+kU8z7VfO4GIXA7cARQaY6JOyOvJRNO3b198Ph+hUOiEJNOtWzfXkkAoFKK5udnVDufw9rgpLS3N9Ztvp6en09Tk3oUd6enpiEi7Wk7Rks3YsWPZunUrlZWVNDY2UlxcTGFhYbyq7p74nd5+F8gTkYEiEgCmA6sjC4jIKOB3WElmTyxBPdlHE/7lDwaDDBo0iI8//piLL77YtUOa5uZmysvLGTFiBN26dXNlHQCVlZVkZmbSt29f19Zx8OBB9uzZw/nnO7lRljPBYJDy8nLGjBnjauJ///33Oe+88zjjjDMcvfeNN95g5cqVPPvsszz//PMnLPP7/SxdupQrr7ySUCjE7NmzGT58eDyrHn9xHOtkjAmKyK3AK1g97o8bYz4SkbuBMmPMaqxDpe7A8/b+3WGMaTMbS5Rft6Scd6KhoYHS0lKam5sJBoOkpaW5ergRTmpunz4P98+49eUEa1tCoZDr1xx1xrY0NzfT3Nzcrm2ZP38+n3/+OcYYevfuTe/evWO6xsYlHfqQJLPAMKwsekGAMnnPGFPQkfW1hycTjVIppmOJpluBYUiMiaY8MYnGk4dOSqkWdFClUspVOnpbKeU6vfGVUsp12qJRSnUKTTRKKVd54MZXmmiU8joPTCDnySEISsXT/v37mTRpEnl5eUyaNIkDBw6cstyKFSvIy8sjLy+PFStWnLS8sLCQ/Px8t6t7sjjf+MoNmmjUaW/RokVMnDiRrVu3MnHixFPeVW///v0sXLiQd955h02bNrFw4cITEtIf/vAHunfv3pnV/lJ47u043PjKLZpo1GmvpKSEwsJCJk2axLJly1i6dOlJrZpXXnmFSZMm8cILL3DRRRexb98+fvrTnwKwZ88eioqKKC0tZdu2bcyff9LdL92n9wxWKrnt3r2bn//85/z1r38lLS2NY8eOndSq2b59Oxs2bKCoqIiePXsyc+ZMnnvuOQ4cOMBVV111fExX37592bhxIy+/HPVeUPFlYnwkiCYadVq4/PLLyc/PP+lRUlKCMYaSkhJeeOEFKioqMMbw7LPPnvD+TZs2cezYMWbPns2PfvQjNmzYQF5eHgsXLqSiooKdO3eyYsUKPvvsM0aOHElNTU2CtjQ56aBKddrr378/O3fu5MUXX+SWW26hqqoKESEY/LL39MILL6SpqQm/3099fT2ffvopV155JTt27KCysvL4/YqCwSB+v58tW7Ycv+9wDDo2qFIKDMQ4qJLEDKrUFo06bbTWqjn77LMxxnDLLbcwdepURo0aRSgUon///scPoZqamqiqqiI9PZ0333wTYwzr169nx44d3HzzzVRWVvL3v/8dn8+Hz+fjwgsv5Ne//nWCtzh56HU06rSxbt26U76+fPlyioqK2LlzJ++99x6VlZUEAgGmTJnC/fffT2FhIT6fj169erFlyxZGjBhBRkYGxhgGDBjAxo0beeCBB7j22mtJS0sjGAzyzW9+sxO3LHzaKXlpi0ad9oYMGcKZZ55Jfn4+v/zlLxERBgwYgN/vZ968eZSUlJCTk0OXLl3IysoiPz+f9PR0QqEQl19+OR988AF33HEHoVCIc845h969ezNixIhO3IL43jTYDZpo1Glv7NixGGPYsWMHV199NYcOHeJ73/seYB0yPfXUUxQWFvLFF18wd+5cNmzYQDAY5Ctf+QpTp06lqamJe+65h9LSUrZv307Pnj3ZvHlzJ25B8l+xp4lGnfb8fj/z58/nyJEjGGP46le/SteuXamoqODAgQNMmDCBOXPmEAgEWLx4MYFAgLVr19Lc3MzXvvY1unfvTk5ODp999hmzZs3i448/ZvTo0Z24BcnfotE+GqWAuXPn8vDDD3P//ffzyCOPUFxczDPPPMPq1dYEABkZGfzsZz/jgw8+oLKyktdff50JEyYQCoUQEbKysujZsycbNmxgwIABHDx4kLS0NDIyMrj11ltdrn3yj6rURKMUX85+MHfuXKqqqvjhD39IXl4eixcv5mc/+xkAc+bM4brrrmPz5s2Ul5dTWlrKqlWryM7O5pNPPiEvL49AIADAd7/7Xc4666xOSDLw5eTbyUuvo1GqhTVr1nDbbbcdn27ljjvu4K677qKgoIDCwkIaGhq47rrrKC8vJysri+Li4pOumVmwYAHdu3dn3rx5sayyg9fR5Bt4PnpBAIbpLAhKnaY6mGiGG1gZY+mv6iwISqn2SP57eWqiUcrzkr8zWE9vK9UJHnnkEUaOHMnIkSMZOHAg48ePj2P05L+ORvtolOpETU1NTJgwgR//+Md8+9vfDr/cwT6a8w38V4ylJ2kfjVKpbu7cuUyYMCEyycRB8h86aaJRqpM88cQTbN++naVLl7oQXTuDlTrtvffee/z6179mw4YNpKXFu2tUWzRKKWDp0qXs37//eCdwQUEBjz76aJyia6JRSmHd88Y9eh2NUsp1yX/jK000SnmeHjoppVyX/IdOemWwUp4X3xtfichkEdkiIttE5KTZ8ESki4g8ay9/R0QGRIupiUYpz4vfEAQR8QEPA98AhgEzRGRYi2JzgAPGmMHAg8CvosXVRKOU58V18u2LgG3GmE+NMY1AMXBVizJXASvsv1cBEyU8VWcrovXRdGgMhlKqM+x6BRb0jrFwhohEzja3zBizLOJ5DlAd8bwGGNcixvEyxpigiBwCzgT2trZS7QxWyuOMMZMTXYdo9NBJKRWpFugX8TzXfu2UZUTED5wB7GsrqCYapVSkd4E8ERkoIgFgOrC6RZnVwPX239OA10yU+83ooZNS6ji7z+VW4BXABzxujPlIRO4Gyowxq4HHgKdEZBuwHysZtSnaja+UUqrD9NBJKeU6TTRKKddpolFKuU4TjVLKdZpolFKu00SjlHKdJhqllOv+P+ZPpYSPrXeGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solving time step:  2\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARoAAAD9CAYAAABjhbg0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAApNElEQVR4nO3df3wU1b3w8c83+4NA+CUGNCYgIFF+BMtP0du++iiIUtuGKtRCwaqQ63189BavD+1DS2vF9ipXawWL91qqIupLonLbBhWxBeuVohEjqUXTImgCSUAQCAjk1+7mPH/MBJdIkp2wk93ZfN+v177I7pw5c2aG/e6ZM+fMEWMMSinlprREF0Aplfo00CilXKeBRinlOg00SinXaaBRSrlOA41SynUaaNogIsdFZGgry24Skb90dpmU8qKEBxoR+a6IlNhf6n0i8oqIfOUM8jMiMqzFZ71E5FciUiEiJ0Rkj4isFZFJbeVljOlpjPm4g+UIisjdIrLT3maFiDwhIoM7kp9b7DI+k+hyuElEBtv/L/yJLktXldBAIyJ3AsuAe4FzgEHAfwLTO5DXaf8TiUg34DVgNPANoDcwAigEvuYkL4fWAvnAd4E+wJeAd4EpcchbxZkGIZcZYxLywvryHQe+3UaaS4C3gCPAPmAFEIxaboDbgJ1AOfCG/dkJO+/vAAX2uhntlOeUvKI+G2b/fTawDvgM2Ar8HPhLK3ldCdQBA9vY3nl2foeBXcA/Ry27G3gBeAY4BmwHLgR+BBwAKoGrotK/Dtxnl+szoAjoZy+7HKhqse0Ku4zTgEYgZB+v96LOzeP2casGfgH4WtmP7sBqoAb4O/DD6O3Z+/nfwKf2Ofp+1LJuWD80e+3XMqBbdLnt/A7YZfkWcA3woX3cfhyVVxqwCPgIOAQ8H3UM9tjn8rj9ugy4CdgCPGSnv9fOc3RUngOAWqB/or4nqfJK3Iat/+RhwN9GmvHApYAfGGz/R74jarkB/gT0A7pHfTYsKk0h8GQM5WkzLzuf54EMIM/+ArYWaJYC/9PO9t7Aqr2lA2PsL+Jke9ndQD1wtb3vT9lf0sVAAPhn7GBop3/dLk+eXb7/Bp6xl11OK4EmalvPtFj+e+A3dl4DsALYv7S1r8BZQA7wt+bt2V/+d4G7gCAwFPgYuNpefg9QbG+jP/Am8POocoftdZv3+VPgWaAXMAormA+x0y+w88rBCmC/AdbYywbb59IfVe6b7Pz/1T7G3e3z8R9RaRYALyb6S5oKr8RtGOYAnzhc5w7g91HvTfOXs8Vn0YFmI7A06v0YrBrSZ8COWPICfFi/+sOjlt1L64Hmt0BhG/sxEIgAvaI+uw87INpf/j9FLfsm1i+xz37fyy5bX/v96y32cSRWTcWHw0CDdQnbgB1s7c9mA39uZV9OBg77fQGfB5pJwJ4W6X8ErLL//gi4JmrZ1UCF/fflWIGk5T5Pikr/LvAt+++/A1OilmXZ56z5R+p0gaZl2SZh1X7Efl8CXJ+o70gqvRJ5XXoIyBQRvzEmfLoEInIh8CtgAtAD6z/Nuy2SVcawnazmN8aYvwJ9ReRK4LEY8+pvbzt6+e52tnlhG8vPAw4bY461yG9C1Pv9UX/XAQeNMZGo9wA9sYImpylbAMhsowytOd9ed5+INH+WRuvH5rwWy6L/Ph84T0SORH3mAzZHrRt9HHfbnzU7dJp9bnlcekZt6/ci0hS1PIIVOFtzyj4ZY94WkVrgchHZh/Ujs66N9VWMEtkY/BbWL+e32kjzX8A/gFxjTG/gx4C0SNPe8PNNwFUikhFDmVrL61OsavbAqM8GtZHPRuASEclpZfleoJ+I9GqRX3UMZWxNy7KFgINY7VU9mheIiA8rcDZruc+VWOcl0xjT1371NsaMamW7+7AuV05XjkqsS7y+Ua9exphr7OV7sQJEdLn3trmXrasEvtZiW+nGmOrT7GOz032+GpgL3ACsNcbUd7A8KkrCAo0x5ijW9fcjIvItEekhIgER+ZqI3G8n64V1iXNcRIYDt8aQ9X6stoBmT2F9GX4vInki4hORdE6tPbRX1gjwO+Buu5wjgRvbSL8Rq73n9yIyXkT89i32/y0i84wxlVjtEfeJSLqIXAzMx2r87ai5IjJSRHpgtX2stcv9IZAuIl8XkQDwE6w2jGb7gcEikmaXfR/wR+BBEektImkicoGI/K9Wtvs88CMROUtEsoHbo5ZtBY6JyP8Tke72sc8TkYn28jXAT0Skv4hkYv1/6OgxeBT4dxE5H8DOs/nu5adAE6f+v2jNM8C1WMHmqQ6WRbWQ0NvbxpgHgTux/vN/ivWrdDvwBzvJQqzbw8ew2j2eiyHbu4HVInJERK63f5GuAMqAl7HbZoCJwPUOins7VjX9E+BJYFU76WcC6+0yHwXexwpuG+3ls7HaDvZiNb7+zA5QHfW0Xa5PsBqYvw8nA/r/wbpMrMaq4VRFrfeC/e8hEdlm//09rMbbMqy7SWuJuvxs4R47v3J739Zi1YiaA/Q3sNrFyrFqWI9h3dUC625WCVYD8nZgm/1ZRyzHusz5o4gcw2oYnmSXoxb4d2CL/f/i0tYysX8EtmHVdja3lk4509zopTxMRF7HatBt2eaUiLLcCswyxrRWA0p6IvIEsNcY85NElyVVaCcldUZEJAvrkuQtIBf4v1j9nTzJ7rl9HTA2wUVJKQkfgqA8L4jVZ+UYVg/sIqz+KJ4jIj/HusR9wBhTnujyJII9TOaAiLzfynIRkYdFZJeI/E1ExsWUr146KaWaichXsfpsPWWMyTvN8muwOjleg9UGttwY0+aYQdAajVIqijHmDayhGK2ZjhWEjDGmGKtPWms3Ck5qr41GqztKua9l3zBHhomY2hjT7oMPsIa3NFtpjFnpYHPZnNrRscr+bF9bK2ljsFIeVwv8S4xp74Z6Y0zMfcjiRQONUh4ndOoXuZpTe3/nEEOPdm2jUcrj0rCGnsfyioN1wPfsu0+XAkft3uRt0hqNUh4nWKNg45KXyBqskfOZIlIF/Kw5e2PMo1i93a/BeoZSLXBzLPlqoFHK4+J56WSMmd3O8uYHxDmigUYpj4tnjcYtGmiU8rhObgzukGQvn1KqHV6o0ehdpxTzzjvvcPHFF1NfX8+JEycYNWoU779/2mErKkV08l2nDtEaTYqZOHEi+fn5/OQnP6Guro65c+eSl/eFISsqhXihRtPeoEodguBBjY2NTJw4kfT0dN588018Pl+ii6TadkZDEHJFzMMxpr0G3tWewSouDh06xPHjxwmFQtTX15OREcvjkpVXaY1GJUR+fj6zZs2ivLycffv2sWKFZ59D1VWcUY3mIhHzmxjTXqE1GhUPTz31FIFAgO9+97tEIhH+6Z/+iddee43JkycnumjKJc2NwclMazRKJd4Z1WhGiphYp44YrzUapVRHaIc9pZTrvNAYrIFGKY/TGo1SynVao1FKuU5I/rtOGmiU8jgBArF+k8NulqR1GmiU8jgR8GugUUq5SQQCST6cTQONUh7nqEaTIElePKVUe0Qg0C3RpWibBhqlvM4DHWlS9gl7GzZs4KKLLmLYsGEsXbo0bvlWVlZyxRVXMHLkSEaNGsXy5cvjlnezSCTC2LFj+cY3vhHXfI8cOcLMmTMZPnw4I0aM4K233opb3g899BCjRo0iLy+P2bNnU19f3/5KbZg3bx4DBgw45aFdhw8fZurUqeTm5jJ16lRqamrOtNipoTnQxPJKkJQMNJFIhNtuu41XXnmFsrIy1qxZQ1lZWVzy9vv9PPjgg5SVlVFcXMwjjzwSt7ybLV++nBEjRsQ1T4AFCxYwbdo0/vGPf/Dee+/FbRvV1dU8/PDDlJSU8P777xOJRCgsLDyjPG+66SY2bNhwymdLly5lypQp7Ny5kylTpsT1B8TzNNB0vq1btzJs2DCGDh1KMBhk1qxZFBUVxSXvrKwsxo0bB0CvXr0YMWIE1dXtzggas6qqKl5++WUKCgrilifA0aNHeeONN5g/fz4AwWCQvn37xi3/cDhMXV0d4XCY2tpazjvvvDPK76tf/Sr9+vU75bOioiJuvPFGAG688Ub+8Ic/nNE2UoYAvhhfCZKSgaa6upqBAz+fHjgnJyeuwaBZRUUFpaWlTJo0KW553nHHHdx///2kpcX31JSXl9O/f39uvvlmxo4dS0FBASdOnIhL3tnZ2SxcuJBBgwaRlZVFnz59uOqqq+KSd7T9+/eTlZUFwLnnnsv+/fvjvg1P0kun1HX8+HFmzJjBsmXL6N27d1zyfOmllxgwYADjx4+PS37RwuEw27Zt49Zbb6W0tJSMjIy4XXrU1NRQVFREeXk5e/fu5cSJEzzzTKxPSOkYEUHkjB7jkjoE6BbjK0E8GWiamppo64Fd2dnZVFZWnnxfVVVFdnZ23LYfCoWYMWMGc+bM4brrrotbvlu2bGHdunUMHjyYWbNm8dprrzF37ty45J2Tk0NOTs7J2tfMmTPZtm1bXPLeuHEjQ4YMoX///gQCAa677jrefPPNuOQd7ZxzzmHfPms++X379jFgwIC4b8OTtEbjjsbGRhobG1sNNhMnTmTnzp2Ul5fT2NhIYWEh+fn5cdm2MYb58+czYsQI7rzzzrjk2ey+++6jqqqKiooKCgsLmTx5ctxqBueeey4DBw5kx44dAGzatImRI0fGJe9BgwZRXFxMbW0txhg2bdrkSmN2fn4+q1evBmD16tVMnz497tvwJA8EmiS/+966PXv2MGjQIILB4Beq0H6/nxUrVnD11VcTiUSYN28eo0aNist2t2zZwtNPP83o0aMZM2YMAPfeey/XXHNNXPJ3069//WvmzJlDY2MjQ4cOZdWqVXHJd9KkScycOZNx48bh9/sZO3Yst9xyyxnlOXv2bF5//XUOHjxITk4OS5YsYdGiRVx//fU8/vjjnH/++Tz//PNxKX9KSPIhCJ58ZnB9fT1vvfUWl156KdXV1VxwwQV6va687Iz+807oJaYkxmY9+R99ZrBjaWlpJ9tiNNioLkt7BneOyspKPvrooy+02UybNs21bXo1b7fz92renuaBu05JHgdjd7qazcGDB13bnlfzdjt/r+btaR6o0SR58ZyprKxk+PA8IpGGk5+1fznlx9nTgD5PH9ulmtP8rXWcXQY63wc38xcJOCxP7Nv4vNxO9zmAMY0Oy+QRGmg6XyTSAOKgDdsI+B2kDwt0d5C+TiDTYZv6QYEsB+vsExjoIH1lB/J3sg8HHR4jsI6T0/Pg9DynquYhCEksJdpolOrS4tyPRkSmicgOEdklIotOs3yQiPxZREpF5G8i0m7fjpSr0SjV5TQ3BscjKxEf8AgwFagC3hGRdcaY6EcU/AR43hjzXyIyElgPDG4rX63RKOV18a3RXALsMsZ8bKxGrUKgZRdsAzQP8OsD7G0vU63RKOV1zhqDM0WkJOr9SmPMyqj32UBl1PsqoOXjCe4G/igi/wpkAFe2t1ENNEqlgti/yQfj0DN4NvCkMeZBEbkMeFpE8owxTWdePKVUcorvXadqYGDU+xz7s2jzgWkAxpi3RCQdyAQOtJapttEo5XXxbaN5B8gVkSEiEgRmAetapNkDTAEQkRFAOvBpW5l6skZz6NAhwuEETbmnVLKJ410nY0xYRG4HXsWqJz1hjPlARO4BSowx64D/C/xWRP4Nq2H4JtPO6GxPBpoePXrQ0NDA7t27E10UpRIvzj2DjTHrsW5ZR392V9TfZcCXneTpyUDTvXt3evTowbFjx6irqyM9PV1HbquuS4cguEdEyMvL4+DBg9TW1tK9e3f7gd4Bh93N/VZ3difp6xymP+g0CPqtbv9O0lc6TO80f0f74PQY2es4PQ9Oz3Oq0kDjvmAwiM/no7a2lm7dugEh98cuneUgfY3DcUhgBY1hDtbZJTDaQfrtHcjf6VgqJ8cIrOPk9Dw4Pc+pLMnHOnk+0AD4fD4yMjKoq6tLdFGU6nxao+k8IkKPHj0SXQyVxCZMmEBmZuYXZsD0vDjedXKL9qNRXUZJSUmrQcbTc317YBYEDTRK4fG5vjXQKOUNnp7r2wNzb6dMG41S8eaZub61MVip1JDUc30L1mijJKaXTkq1wjNzfXvg0kkDjVKt8Mxc39oYrJQ3zJ49m8suu4wdO3aQk5PD448/zqJFi/jTn/5Ebm4uGzduZNGiLzynO3kkeaBJwTaagPtjl2rcHIdkr7PL4Ta2O0zvNH+nY6kcHSN7Hafnwel5bsOaNWtO+/mmTZscbCNBPDDdSgoGmpD7Y5eGOEhfLjDW4bifUoErHKzzZ4FvOkj/Ygfyd7IPpQ6PEVjHyel5cHqeU5XedVJKuc4DQxA00CjldVqjUUq5TgONUsp1GmiUUp1C7zoppVylNRqllOv0rpNSynVao3FHJBJJdBGUSh4aaNxx+PBhjh8/TkVFBcaY5B2+r1Rn0EDjjv79+5ORkQHAiRMn8Pv9BINBe14npboeo3ed3CEiDB48mL179xIKhairq7NrNgH3B0mWO0xf6rDGJX5rfJGT9C86TO8kf8f74PQY2es4PQ9Oz3OKMmnQmOQPvkqJox8IBAgEAnbbTQgyHQy2O9iBydGcDjB0MuARrKDxbw7WeUjg3x2kX9yB/J0O2uzIQFKn58HheU7V6VaMQNgXa22+ydWytCYlAk0zny/J648qoUpKSmJK99BDD/HYY48hIowePZpVq1aRnp68VQYjQsQf61e50dWytEYbNZSKUl1dzcMPP0xJSQnvv/8+kUiEwsLCRBerXRGfL6ZXoqRUjUapeAiHw9TV1REIBKitreW8885LdJHaZBAiST4GQWs0SkXJzs5m4cKFDBo0iKysLPr06cNVV12V6GK1ySCE8cX0ShQNNEpFqampoaioiPLycvbu3cuJEyd45plnEl2sNhmERrrF9EoUDTRKRdm4cSNDhgyhf//+BAIBrrvuOt58881EF6tNzZdOsbwSRQONUlEGDRpEcXExtbW1GGPYtGkTI0aMSHSx2hXPQCMi00Rkh4jsEpHTTv0gIteLSJmIfCAiz7aXpzYGKxVl0qRJzJw5k3HjxuH3+xk7diy33HJLoovVpuY2mngQER/wCDAVqALeEZF1xpiyqDS5wI+ALxtjakSk3Zn1NNAo1cKSJUtYsmRJoosRM+vSKW5f5UuAXcaYjwFEpBCYDpRFpfln4BFjTA2AMeZAe5lqoFHK46zG4GCsyTNFJLrn4kpjzMqo99lAZdT7KmBSizwuBBCRLVjP9rvbGNNmd+sUDDQBa1hBzDowOZqTcT9OxyEBpPmtbv9O0i92mN5J/o73oQPjuzpyHpye5xRlwMml00FjzIQz3KQfyAUuB3KAN0RktDHmSFsrpJgQZDkYA7NPYJiD9Ls6MPmak3FFYAWBIgfrTBf6NOyLOfnRblmO83c8NsrJMQLrODk9D07Pc8qK66VTNTAw6n2O/Vm0KuBtY0wIKBeRD7ECzzutZap3nZTyuDjf3n4HyBWRISISBGYB61qk+QNWbQYRycS6lPq4rUxTsEajVNcTrz4yxpiwiNwOvIrV/vKEMeYDEbkHKDHGrLOXXSUiZUAE+IEx5lBb+WqgUcrj4j3WyRizHljf4rO7ov42wJ32KyYaaJTyOIPQkOTTIGigUcrjvDB6WwONUh6ngUYp1SkS+QiIWGigUcrj4jwEwRXJXTqlVLv00sklNTU11NXVsWvXLhobGxER0tLSdCI51SVZd51iHuuUEJ4MNBkZGQQCAXr16gVAU1MTkUiEpqYmIOCwu7nf6s4eK6dzIjkdVwTg81vd/mPl91vDCtzKvyNjoxzNGwWOzwN+x+c5Zadb0UsndwSDQfx+P+eccw4fffRRi6Uh5/MDjXaQfnsH5jhyMucSwGLnY5f2mT4xp8+So87HRjmdN6ojc1k5PQ8Oz3Os060cOXKEgoIC3n//fUSEJ554gssuuyz2bSWAXjop5TELFixg2rRprF27lsbGRmpraxNdpDZpG41SHnP06FHeeOMNnnzyScCqPQeDyd3+4YVAo6O3lYpSXl5O//79ufnmmxk7diwFBQWcOHEipnXvuusuli1bdvL94sWLWb58uUsl/VzzEIRYXomigUapKOFwmG3btnHrrbdSWlpKRkYGS5cujWndefPm8dRTTwHWDYrCwkLmzp3rZnEBnQVBKc/JyckhJyeHSZOsp1fOnDmTbdu2xbTu4MGDOfvssyktLeWPf/wjY8eO5eyzz3azuCcle6DRNhqlopx77rkMHDiQHTt2cNFFF7Fp0yZGjhwZ8/oFBQU8+eSTfPLJJ8ybN8/Fkn4unrMguEUDjVIt/PrXv2bOnDk0NjYydOhQVq1aFfO61157LXfddRehUIhnn213uqO40H40SnnQmDFjYu5z01IwGOSKK66gb9+++HydV8tI9rtOGmiUiqOmpiaKi4t54YUXOm2bDqdbSQhtDFYqTsrKyhg2bBhTpkwhNze307bb3EYTyytRUq9GIwHn8wNtd3GOI6dzLoHjsUs+vzWsIPb8fc7GRjndh47MZeX0PHRkHiiXjRw5ko8/bnMyAFdoG00imK43r1NE53VqX0rP66RtNEopl3lhCIIGGqU8TvvRKKVcZ9110ulWlFIu0ksnpVSn0ECjlHKVttEopVyn/WiUUq7zwhAEDTRKeZxeOnWSpqYmwuEwkUgk0UVRKiH00skFoVCIUCjEBx98wPHjx0lLS8Pn81kPkRad1ynu+XfBeZ0ikQgTJkwgOzubl156ycF2Op/e3nbJsWPHaGpqIisri6NHWwwmNCHIdDAG5qDz+YEY6yB9aQfnOHI6tsjpvEtO83c6l5WTYwTWcXJ6HhyeZyfPmFm+fDkjRozgs88+i30bCeKFQOPJx0T069ePbt260a9fv0QXRaWgqqoqXn75ZQoKChJdlJjpYyKU8pg77riD+++/n2PHjiW6KDFpIi3phyB4skajlFteeuklBgwYwPjx4xNdFEfiOQuCiEwTkR0isktEFrWRboaIGBGZ0F6eWqNRKsqWLVtYt24d69evp76+ns8++4y5c+fyzDPPJLporYpnG42I+IBHgKlAFfCOiKwzxpS1SNcLWAC8HUu+WqNRKsp9991HVVUVFRUVFBYWMnny5KQOMgCGuLbRXALsMsZ8bIxpBAqB6adJ93PgP4D6WDLVQKOU51lDEGJ5AZkiUhL1uqVFZtlAZdT7Kvuzz7cmMg4YaIx5OdYS6qWTUq24/PLLufzyyxNdjHY5vHQ6aIxpt02lNSKSBvwKuMnJehpolPI4g9AQv7FO1cDAqPc59mfNegF5wOsiAnAusE5E8o0xrXZU0kCjlMfFefT2O0CuiAzBCjCzgO+e3JYxR4HM5vci8jqwsK0gAxpolEoJ8brrZIwJi8jtwKuAD3jCGPOBiNwDlBhj1nUkXw00SnlcvIcgGGPWA+tbfHZXK2kvjyXPFAw0AWv8Usw6MBFZqcuTqTkdxOh0greODJJ0tA8Oj1HzOk7Pg9PznKIMQqQpucc6peDRD0F3B4Pt6gTOcpC+RmCIg/TlHRxg6HSSOqeDHp3m73QgqZNjBNZxcnoenJ7nFGWahIb65B6CkIKBRqmuxRghEtYajVLKTQYNNEopdxkjhEMaaJRSrhKaIsn9VU7u0iml2mcAvXRSSrmqSaA+ub/KyV06pVRswokuQNtSKtDodCuqS7IeSJPUUuJ5NMYY6uvrqa+P6Rk8qouaMGEC06ZNS3Qx4q850MTyShBP12iMMYRCIRoaGujWrRvduiV370iVWLFMt1JZWcn3vvc99u/fj4hwyy23sGDBgk4o3RkwQCjRhWibZwNNU1MTJSUlhMNhMjIysJ+NAQQcdjf3W93ZnaQvd3ncj9MJ2JyORXI8wZvTfXB6jOx1nJ4Hp+c5llR+Pw8++CDjxo3j2LFjjB8/nqlTpzJy5EgH2+pkBmhIdCHa5slAc+DAAerq6sjLy6OsrKzF0hD4HYyBCXdgzIzTMTlOJkYDa3DhMAfr7BIY7SD99g7k73RyNyfHCDo2dsnpeY5BVlYWWVnWrJ+9evVixIgRVFdXJ3+gSfI2Gk8Gmn79+pGRkUHfvn0TXRSVwioqKigtLWXSpEmJLkrbNNC4w+/3ZLGVhxw/fpwZM2awbNkyevfunejitE0DjVLeEwqFmDFjBnPmzOG6665LdHHap4FGKW8xxjB//nxGjBjBnXfemejixC7JA01K9KNRKl62bNnC008/zWuvvcaYMWMYM2YM69evb3/FRGrCmsYtlleCaI1GqShf+cpXMMbhHbNE00snpZTrNNAopVyngUYp1Sk00CilXKU1mkQIxNzd3NKBMTNOx+Q4mq/IXmeXw21sd5jeaf5O51xydIzsdZyeB6fnOVU1AXWJLkTbUvDoh0Ac3DUwHRgz43RMTqbDuxgHBbIcrLOvA2ORnObvZB8OOjxG0LGxS07Pc6oyQJI/iikFA41SXZBeOimlXKVtNEop12mgUUq5rnkIQhLTQKNUKtAajVLKVXrppJRynQceTp4yj4kwxlBbW5voYqgkltLTrURifCVISgSa5iCjj/hUbSkpKWHDhg3tptuwYQMXXXQRw4YNY+nSpZ1QsjMU53mdRGSaiOwQkV0isug0y+8UkTIR+ZuIbBKR89vL0/OBpjnIBINBAoFAooujPC4SiXDbbbfxyiuvUFZWxpo1a04z00aSMVhDEGJ5tUNEfMAjwNeAkcBsEWk5BUQpMMEYczGwFri/vXw9XQVobGyktraWbt264ff7McbQs2dfjh930t3c6dgop/NGBawu+Y4ErG7/sZKAs7FI4jB/x/vg9BjZ6zg9D46GFQQ4duwYvXr1ajPV1q1bGTZsGEOHDgVg1qxZFBUVJf90K/G7LLoE2GWM+RhARAqB6cDJaGuM+XNU+mJgbnuZejbQNE8gFx1kfD4fa9cWEgqF6N69e9SkcvETCoWIRCKkp6fHPW+wgqeIuFY7q6+vJxAI4PP5XMm/traW9PR00tLiX1lurr12797dUf6LFi3i0KFDDBs2jOzsbAYMGNDqJVR1dTUDBw48+T4nJ4e33377jMvuKmd3nTJFJHrKzpXGmJVR77OByqj3VUBb883MB15pb6OevHSqq6ujtraW4cOHnxJkGhoaXA0yxhgaGxtdnXo3HA67FgQARMTVR1UGg0EaGxtdyVtESE9Pp66uztE+LF26lN/+9rdkZWVRXV3NgQMHUqtR2FkbzUFjzISo18rT5hkDEZkLTAAeaC+tJwPNsWPH6N69O/369TslyDQ0NLgWZAAaGhoIBoOu5W+MwRjjSm2gmYjQ1NTkWv5+v59IJOJaMPP5fASDQerrnXeF/dWvftVusMnOzqay8vMf9KqqKrKzs8+ozK5rvr0dy6t91cDAqPc59menEJErgcVAvjGm3Ql5PRloBgwYgM/nIxKJnBJkevTo4VoQiEQiNDU1udrg3Lw/bkpLS3P94duBQIBQyL2OHYFAABHpUM2pvWAzceJEdu7cSXl5OY2NjRQWFpKfnx+vorsnfre33wFyRWSIiASBWcC66AQiMhb4DVaQORBLpp5so2n+5Q+HwwwdOpQPP/yQSy+91LVLmqamJkpLSxk9ejQ9evRwZRsA5eXlZGRkMGDAANe2ceTIEQ4cOMCFF17o2jbC4TClpaWMHz/e1cD/3nvvccEFF9CnTx9H677xxhusWbOG5557jhdeeOGUZX6/nxUrVnD11VcTiUSYN28eo0aNimfR4y+OY52MMWERuR14FfABTxhjPhCRe4ASY8w6rEulnsAL9vndY4xpMxpLO79uSTnvRH19PcXFxTQ1NREOh0lLS3P1cqM5qLl9+7y5fcatLydY+xKJRFzvc9QZ+9LU1ERTU1OH9mXRokV8+umnGGPIzMwkMzMzpj42LjmjgyQZEwwjS9pPCFAi7xpjJpzJ9jrCk4FGqRRzZoGmxwTD8BgDTWliAo0nL52UUi3ooEqllKt09LZSynX64CullOu0RqOU6hQaaJRSrvLAg6800CjldR6YQM6TQxCUiqfDhw8zdepUcnNzmTp1KjU1NadNt3r1anJzc8nNzWX16tVfWJ6fn09eXp7bxf2iOD/4yg0aaFSXt3TpUqZMmcLOnTuZMmXKaZ+qd/jwYZYsWcLbb7/N1q1bWbJkySkB6Xe/+x09e/bszGJ/rnnu7Tg8+MotGmhUl1dUVER+fj5Tp05l5cqVrFix4gu1mldffZWpU6fy4osvcskll3Do0CF+/OMfA3DgwAEKCgooLi5m165dLFr0hadfuk+fGaxUctu/fz8/+9nP+Otf/0paWhoNDQ1fqNXs3r2bzZs3U1BQQO/evZkzZw7PP/88NTU1TJ8+/eSYrgEDBrBlyxZeeaXdZ0HFl4nxlSAaaFSXcOWVV5KXl/eFV1FREcYYioqKePHFFykrK8MYw3PPPXfK+lu3bqWhoYF58+bxgx/8gM2bN5Obm8uSJUsoKytj7969rF69mk8++YQxY8ZQVVWVoD1NTjqoUnV5gwYNYu/evbz00kvcdtttVFRUICKEw5+3nl588cWEQiH8fj91dXV8/PHHXH311ezZs4fy8vKTzysKh8P4/X527Nhx8rnDMTizQZUywUCMgypJzKBKrdGoLqO1Ws25556LMYbbbruNGTNmMHbsWCKRCIMGDTp5CRUKhaioqCAQCPCXv/wFYwybNm1iz5493HrrrZSXl/P3v/8dn8+Hz+fj4osv5pe//GWC9zh5aD8a1WVs3LjxtJ+vWrWKgoIC9u7dy7vvvkt5eTnBYJBrrrmGBx54gPz8fHw+H3379mXHjh2MHj2a9PR0jDEMHjyYLVu28OCDD3L99deTlpZGOBzm61//eifuWfNtp+SlNRrV5Q0fPpyzzz6bvLw8fvGLXyAiDB48GL/fz8KFCykqKiI7O5tu3brRr18/8vLyCAQCRCIRrrzySrZv387ixYuJRCKcd955ZGZmMnr06E7cg/g+NNgNGmhUlzdx4kSMMezZs4drr72Wo0eP8p3vfAewLpmefvpp8vPz+eyzz1iwYAGbN28mHA5z1llnMWPGDEKhEPfeey/FxcXs3r2b3r17s23btk7cg+TvsaeBRnV5fr+fRYsWcfz4cYwxfOlLX6J79+6UlZVRU1PD5MmTmT9/PsFgkGXLlhEMBtmwYQNNTU18+ctfpmfPnmRnZ/PJJ58wd+5cPvzwQ8aNG9eJe5D8NRpto1EKWLBgAY888ggPPPAAjz76KIWFhTz77LOsW2dNAJCens5Pf/pTtm/fTnl5Oa+//jqTJ08mEokgIvTr14/evXuzefNmBg8ezJEjR0hLSyM9PZ3bb7/d5dIn/6hKDTRK8fnsBwsWLKCiooLvf//75ObmsmzZMn76058CMH/+fG644Qa2bdtGaWkpxcXFrF27lqysLD766CNyc3MJBoMAfPvb3+acc87phCADn0++nby0H41SLaxfv5477rjj5HQrixcv5q677mLChAnk5+dTX1/PDTfcQGlpKf369aOwsPALfWbuvvtuevbsycKFC2PZ5Bn2o8kz8EL7CQEYqbMgKNVFnWGgGWVgTYypv6SzICilOiL5n+WpgUYpz0v+xmC9va1UJ3j00UcZM2YMY8aMYciQIVxxxRVxzD35+9FoG41SnSgUCjF58mR++MMf8s1vfrP54zNso7nQwH/GmHqqttEoleoWLFjA5MmTo4NMHCT/pZMGGqU6yZNPPsnu3btZsWKFC7lrY7BSXd67777LL3/5SzZv3kxaWrybRrVGo5QCVqxYweHDh082Ak+YMIHHHnssTrlroFFKYT3zxj3aj0Yp5brkf/CVBhqlPE8vnZRSrkv+SyftGayU58X3wVciMk1EdojILhH5wmx4ItJNRJ6zl78tIoPby1MDjVKeF78hCCLiAx4BvgaMBGaLyMgWyeYDNcaYYcBDwH+0l68GGqU8L66Tb18C7DLGfGyMaQQKgekt0kwHVtt/rwWmSPNUna1or43mjMZgKKU6w75X4e7MGBOni0j0bHMrjTEro95nA5VR76uASS3yOJnGGBMWkaPA2cDB1jaqjcFKeZwxZlqiy9AevXRSSkWrBgZGvc+xPzttGhHxA32AQ21lqoFGKRXtHSBXRIaISBCYBaxrkWYdcKP990zgNdPO82b00kkpdZLd5nI78CrgA54wxnwgIvcAJcaYdcDjwNMisgs4jBWM2tTeg6+UUuqM6aWTUsp1GmiUUq7TQKOUcp0GGqWU6zTQKKVcp4FGKeU6DTRKKdf9f3NaTtv+6WtwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solving time step:  3\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARoAAAD9CAYAAABjhbg0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAApjElEQVR4nO3df3gU1b348fdnd7MEAkgxoDEBAYnyI1h+BNG2T6+CCKU2VqEWClaFXO/1iy1+lfZLi1WxrXL9UUHxXktVRHwkKrc1qIgteLlSFDGQWjEtgiaSBATlNyQh2d3z/WN2cQkk2Ul2sjvL5/U888DunDlzZif72TNnzpkjxhiUUspJnkQXQCmV+jTQKKUcp4FGKeU4DTRKKcdpoFFKOU4DjVLKcRpomiEiR0WkXxPrbhKRv7Z3mZRyo4QHGhH5kYiUhL/Uu0XkDRH5VhvyMyLSv9F7XUTkdyJSISLHRGSniKwQkVHN5WWM6WyM+bSV5fCLyL0isj28zwoReUZE+rQmP6eEy/h8osvhJBHpE/678CW6LGeqhAYaEbkDWADcD5wD9Ab+E7imFXmd9o9IRDoAbwFDgKuBrsBAoAj4jp28bFoBFAA/As4Cvg5sBsbEIW8VZxqEHGaMSciC9eU7CvygmTSXAO8CB4HdwCLAH7XeADOB7UA58Hb4vWPhvH8IFIa3zWihPCflFfVe//D/zwZWAoeBTcCvgb82kdeVQC3Qq5n9nRfObz+wA/jXqHX3Ai8DzwNHgA+BC4FfAHuBSuCqqPTrgAfC5ToMFAPdw+suB6oa7bsiXMbxQD3QEP68Pog6N0+HP7dq4DeAt4nj6AgsBQ4A/wB+Hr2/8HH+N/BF+Bz9NGpdB6wfml3hZQHQIbrc4fz2hsvyfWAC8HH4c/tlVF4eYA7wCbAPeCnqM9gZPpdHw8tlwE3ABuDRcPr7w3kOicqzJ1AD9EjU9yRVlsTt2PojDwC+ZtKMAC4FfECf8B/y7VHrDfAXoDvQMeq9/lFpioBnYyhPs3mF83kJyADywl/ApgLNfOB/W9jf21i1t3RgaPiLODq87l6gDhgXPvbnwl/SuUAa8K+Eg2E4/bpwefLC5ftv4PnwustpItBE7ev5Ruv/BPw+nFdPrAD2b80dK/A1IAf4e2R/4S//ZuBuwA/0Az4FxoXX3wdsDO+jB/AO8OuocgfC20aO+QvgBaALMBgrmPcNp58VzisHK4D9HlgeXtcnfC59UeW+KZz/T8Kfccfw+fiPqDSzgFcT/SVNhSVxO4apwOc2t7kd+FPUaxP5cjZ6LzrQrAHmR70eilVDOgxsiyUvwIv1qz8gat39NB1o/gAUNXMcvYAg0CXqvQcIB8Twl/8vUeu+h/VL7A2/7hIuW7fw63WNjnEQVk3Fi81Ag3UJe5xwsA2/NwX4nyaO5UTgCL8u5KtAMwrY2Sj9L4Al4f9/AkyIWjcOqAj//3KsQNL4mEdFpd8MfD/8/38AY6LWZYXPWeRH6nSBpnHZRmHVfiT8ugS4PlHfkVRaEnldug/IFBGfMSZwugQiciHwOyAf6IT1R7O5UbLKGPaTFXlhjPkb0E1ErgSeijGvHuF9R6//rIV9XtjM+vOA/caYI43yy496vSfq/7XAl8aYYNRrgM5YQZPTlC0NyGymDE05P7ztbhGJvOeh6c/mvEbrov9/PnCeiByMes8LrI/aNvpz/Cz8XsS+0xxz48+lc9S+/iQioaj1QazA2ZSTjskY856I1ACXi8hurB+Zlc1sr2KUyMbgd7F+Ob/fTJr/Av4J5BpjugK/BKRRmpaGn68FrhKRjBjK1FReX2BVs3tFvde7mXzWAJeISE4T63cB3UWkS6P8qmMoY1Mal60B+BKrvapTZIWIeLECZ0TjY67EOi+Zxphu4aWrMWZwE/vdjXW5crpyVGJd4nWLWroYYyaE1+/CChDR5d7V7FE2rRL4TqN9pRtjqk9zjBGne38pMA24AVhhjKlrZXlUlIQFGmPMIazr7ydE5Psi0klE0kTkOyLyYDhZF6xLnKMiMgC4NYas92C1BUQ8h/Vl+JOI5ImIV0TSObn20FJZg8AfgXvD5RwE3NhM+jVY7T1/EpERIuIL32L/dxGZboypxGqPeEBE0kXkYmAGVuNva00TkUEi0gmr7WNFuNwfA+ki8l0RSQPuwmrDiNgD9BERT7jsu4E/A4+ISFcR8YjIBSLyL03s9yXgFyLyNRHJBm6LWrcJOCIi/09EOoY/+zwRGRlevxy4S0R6iEgm1t9Daz+DJ4Hfisj5AOE8I3cvvwBCnPx30ZTngWuxgs1zrSyLaiSht7eNMY8Ad2D98X+B9at0G/BKOMlsrNvDR7DaPV6MIdt7gaUiclBErg//Il0BlAGvE26bAUYC19so7m1Y1fTPgWeBJS2knwSsCpf5ELAVK7itCa+fgtV2sAur8fWecIBqrWXhcn2O1cD8UzgR0P8P1mViNVYNpypqu5fD/+4TkS3h//8Yq/G2DOtu0gqiLj8buS+cX3n42FZg1YgiAfpqrHaxcqwa1lNYd7XAuptVgtWA/CGwJfxeayzEusz5s4gcwWoYHhUuRw3wW2BD+O/i0qYyCf8IbMGq7axvKp2yJ9LopVxMRNZhNeg2bnNKRFluBSYbY5qqASU9EXkG2GWMuSvRZUkV2klJtYmIZGFdkrwL5AJ3YvV3cqVwz+3rgGEJLkpKSfgQBOV6fqw+K0ewemAXY/VHcR0R+TXWJe5DxpjyRJcnEcLDZPaKyNYm1ouIPCYiO0Tk7yIyPKZ89dJJKRUhIt/G6rP1nDEm7zTrJ2B1cpyA1Qa20BjT7JhB0BqNUiqKMeZtrKEYTbkGKwgZY8xGrD5pTd0oOKGlNhqt7ijlvMZ9w2zpL2JqYky7Gz7CGt4SsdgYs9jG7rI5uaNjVfi93c1tpI3BSrlcDfBvMaa9F+qMMTH3IYsXDTRKuZzQrl/kak7u/Z1DDD3atY1GKZfzYA09j2WJg5XAj8N3ny4FDoV7kzdLazRKuZxgjYKNS14iy7FGzmeKSBVwTyR7Y8yTWL3dJ2A9Q6kGuDmWfDXQKOVy8bx0MsZMaWF95AFxtmigUcrl4lmjcYoGGqVcrp0bg1sl2cunlGqBG2o0etcpxbz//vtcfPHF1NXVcezYMQYPHszWracdtqJSRDvfdWoVrdGkmJEjR1JQUMBdd91FbW0t06ZNIy/vlCErKoW4oUbT0qBKHYLgQvX19YwcOZL09HTeeecdvF5vooukmtemIQi5IuaxGNNOgM3aM1jFxb59+zh69CgNDQ3U1dWRkRHL45KVW2mNRiVEQUEBkydPpry8nN27d7NokWufQ3WmaFON5iIR8/sY016hNRoVD8899xxpaWn86Ec/IhgM8o1vfIO33nqL0aNHJ7poyiGRxuBkpjUapRKvTTWaQSIm1qkjRmiNRinVGtphTynlODc0BmugUcrltEajlHKc1miUUo4Tkv+ukwYapVxOgLRYv8kBJ0vSNA00SrmcCPg00CilnCQCaUk+nE0DjVIuZ6tGkyBJXjylVEtEIK1DokvRPA00SrmdCzrSpOwT9lavXs1FF11E//79mT9/ftzyrays5IorrmDQoEEMHjyYhQsXxi3viGAwyLBhw7j66qvjmu/BgweZNGkSAwYMYODAgbz77rtxy/vRRx9l8ODB5OXlMWXKFOrq6lreqBnTp0+nZ8+eJz20a//+/YwdO5bc3FzGjh3LgQMH2lrs1BAJNLEsCZKSgSYYDDJz5kzeeOMNysrKWL58OWVlZXHJ2+fz8cgjj1BWVsbGjRt54okn4pZ3xMKFCxk4cGBc8wSYNWsW48eP55///CcffPBB3PZRXV3NY489RklJCVu3biUYDFJUVNSmPG+66SZWr1590nvz589nzJgxbN++nTFjxsT1B8T1NNC0v02bNtG/f3/69euH3+9n8uTJFBcXxyXvrKwshg8fDkCXLl0YOHAg1dUtzggas6qqKl5//XUKCwvjlifAoUOHePvtt5kxYwYAfr+fbt26xS3/QCBAbW0tgUCAmpoazjvvvDbl9+1vf5vu3buf9F5xcTE33ngjADfeeCOvvPJKm/aRMgTwxrgkSEoGmurqanr1+mp64JycnLgGg4iKigpKS0sZNWpU3PK8/fbbefDBB/F44ntqysvL6dGjBzfffDPDhg2jsLCQY8eOxSXv7OxsZs+eTe/evcnKyuKss87iqquuikve0fbs2UNWVhYA5557Lnv27In7PlxJL51S19GjR5k4cSILFiyga9euccnztddeo2fPnowYMSIu+UULBAJs2bKFW2+9ldLSUjIyMuJ26XHgwAGKi4spLy9n165dHDt2jOefj/UJKa0jIoi06TEuqUOADjEuCeLKQBMKhWjugV3Z2dlUVlaeeF1VVUV2dnbc9t/Q0MDEiROZOnUq1113Xdzy3bBhAytXrqRPnz5MnjyZt956i2nTpsUl75ycHHJyck7UviZNmsSWLVvikveaNWvo27cvPXr0IC0tjeuuu4533nknLnlHO+ecc9i925pPfvfu3fTs2TPu+3AlrdE4o76+nvr6+iaDzciRI9m+fTvl5eXU19dTVFREQUFBXPZtjGHGjBkMHDiQO+64Iy55RjzwwANUVVVRUVFBUVERo0ePjlvN4Nxzz6VXr15s27YNgLVr1zJo0KC45N27d282btxITU0NxhjWrl3rSGN2QUEBS5cuBWDp0qVcc801cd+HK7kg0CT53fem7dy5k969e+P3+0+pQvt8PhYtWsS4ceMIBoNMnz6dwYMHx2W/GzZsYNmyZQwZMoShQ4cCcP/99zNhwoS45O+kxx9/nKlTp1JfX0+/fv1YsmRJXPIdNWoUkyZNYvjw4fh8PoYNG8Ytt9zSpjynTJnCunXr+PLLL8nJyWHevHnMmTOH66+/nqeffprzzz+fl156KS7lTwlJPgTBlc8Mrqur49133+XSSy+lurqaCy64QK/XlZu16Y83v4uYkhib9eR/9ZnBtnk8nhNtMRps1BlLewa3j8rKSj755JNT2mzGjx/v2D7dmrfT+bs1b1dzwV2nJI+DsTtdzebLL790bH9uzdvp/N2at6u5oEaT5MWzp7KykgGD8gg2HD/xXouXU+IDY+NpQFHpY7pU8/ggZPNpQx6fvctAu/twOH/xpNn7TCHm83Ci3HbPmycNE6y3Vya30EDT/oINx2GKjTbs5QI32Ei/TKDQRvqnBH5is039cYHZNrZ5WOAXNtI/0Ir87RzD4zY/I7A+J7vnwe55TlWRIQhJLCXaaJQ6o8W5H42IjBeRbSKyQ0TmnGZ9bxH5HxEpFZG/i0iLfTtSrkaj1Bkn0hgcj6xEvMATwFigCnhfRFYaY6IfUXAX8JIx5r9EZBCwCujTXL5ao1HK7eJbo7kE2GGM+dQYUw8UAY27YBsgMsDvLGBXS5lqjUYpt7PXGJwpIiVRrxcbYxZHvc4GKqNeVwGNH09wL/BnEfkJkAFc2dJONdAolQpi/yZ/GYeewVOAZ40xj4jIZcAyEckzxoTaXjylVHKK712naqBX1Ouc8HvRZgDjAYwx74pIOpAJ7G0qU22jUcrt4ttG8z6QKyJ9RcQPTAZWNkqzExgDICIDgXTgi+YydWWNZt++fQQCCZpyT6lkE8e7TsaYgIjcBryJVU96xhjzkYjcB5QYY1YCdwJ/EJH/i9UwfJNpYXS2KwNNp06dOH78OJ999lmii6JU4sW5Z7AxZhXWLevo9+6O+n8Z8E07eboy0HTs2JFOnTpx5MgRamtrSU9P15Hb6szlgiEIrn4ezTe+8Q3WrVtHQ0MDHTt2xOPxMPrKcRBqiD2zNox1ciQ9tGrskqPp2+OY2+E8GDt/F+2rbc+jyRFTMjPGHf1Sn0fTKn6/H6/XS01NDR06dLCCjN0xMzfbSL+kFeN+7IxDAmss0jwb29wj8Fsb6ee2In+7Y6laM77L7nmwe55TWZKPdXJ9oAHwer1kZGRQW1ub6KIo1f5ccOmU5MWLnYjQqVOnRBdDJbH8/HwyMzNPmQHT9eJ418kp2o9GnTFKSkqaDDKunuvbBbMgaKBRCpfP9a2BRil3cPVc3y6Yeztl2miUijfXzPWtjcFKpYaknutbsEYbJTG9dFKqCa6Z69sFl04aaJRqgmvm+tbGYKXcYcqUKVx22WVs27aNnJwcnn76aebMmcNf/vIXcnNzWbNmDXPmnPKc7uSR5IHG9WOd3nnnnZPWJd1Yp1bO65RUY52cHhvVmm10rNMJ+ReKKVkU447G6Vin+Ag1ODvvkt2xSw8IPGIzXt8p8LSNbWYIFNlIP7kV+ds5hjtbOb7LybmjnkrShtx40LtOSinHuWAIggYapdxOazRKKcdpoFFKOU4DjVKqXeiDr5RSjtIajVLKcXrXSSnlOK3ROCMYDCa6CEolDw00zti/fz9Hjx6loqICY0zyDt9Xqj1ooHFGjx49yMjIAODYsWP4fD78fj8ej44RVWcmk+R3nVJiUGVDQwP19fWICFcXXOvsoEqnByS2ZhuvD4JJdAytOWYdVNlqI4aLeeevsaVNz9BBla2WlpZGWlqa1XYTakiuQZJ32hzwCNagx802thkhUF0Xe/rsdPv52x202ZqBpE5OUve4pOx0K0Yg4I21Nh9ytCxNSYlAE+H1Jnn9USVUSUlJTOkeffRRnnrqKUSEIUOGsGTJEtLTk/dZmUaEoC/Wr3K9o2VpijZqKBWlurqaxx57jJKSErZu3UowGKSoqCjRxWpR0OuNaUmUlKrRKBUPgUCA2tpa0tLSqKmp4bzzzkt0kZplEIJJPgZBazRKRcnOzmb27Nn07t2brKwszjrrLK666qpEF6tZBiGAN6YlUTTQKBXlwIEDFBcXU15ezq5duzh27BjPP/98oovVLINQT4eYlkTRQKNUlDVr1tC3b1969OhBWloa11133SnPpU42kUunWJZE0UCjVJTevXuzceNGampqMMawdu1aBg4cmOhitSiegUZExovINhHZISKnnfpBRK4XkTIR+UhEXmgpT20MVirKqFGjmDRpEsOHD8fn8zFs2DBuueWWRBerWZE2mngQES/wBDAWqALeF5GVxpiyqDS5wC+AbxpjDohIizPraaBRqpF58+Yxb968RBcjZtalU9y+ypcAO4wxnwKISBFwDVAWleZfgSeMMQcAjDF7W8pUA41SLmc1BvtjTZ4pItE9FxcbYxZHvc4GKqNeVwGjGuVxIYCIbMB6tt+9xphmu1unXqDxpFnDCmJO77O6s9tJf6eN9F6f1SXfDq/P6vYfK5/PGlbgVP52j8HuZxTZxs55EJ/985yiDNi5dPoyDmOdfEAucDmQA7wtIkOMMQeb2yC1hBpgto0xMA8LzLOR/p5WTL5mZ1wRwAjhArM15uSfSB43nPSj1Lxlcovt/G2PjbLzGYH1Odk9D3bPc8qK66VTNdAr6nVO+L1oVcB7xpgGoFxEPsYKPO83lanedVLK5eJ8e/t9IFdE+oqIH5gMrGyU5hWs2gwikol1KfVpc5mmXo1GqTNQvPrIGGMCInIb8CZW+8szxpiPROQ+oMQYszK87ioRKQOCwM+MMfuay1cDjVIuF++xTsaYVcCqRu/dHfV/A9wRXmKigUYplzMIx5N8GgQNNEq5nBtGb2ugUcrlNNAopdpFIh8BEQsNNEq5XJyHIDgiuUunlGqRXjo55MCBA9TW1rJjx44T06x4PB6dSE6dkay7TjGPdUoIVwaajIwM0tLS6NKlCwChUIhgMEgoFLLGOtnpbu7xWd3Z7aSfYXOckJ1xRQA+r9XtP0bi87BMbDzKwGb+to/B7mcU2cbuebB5nlN2uhW9dHKG3+/H5/Nxzjnn8Mknn5y8MtRgf36g39pIP7cVcxzZmXMJIDvd9tgl07iTeDOkIGg7f9vzRrVmLiu758HmeY51upWDBw9SWFjI1q1bERGeeeYZLrvsstj3lQB66aSUy8yaNYvx48ezYsUK6uvrqampSXSRmqVtNEq5zKFDh3j77bd59tlnAav27Pcnd/uHGwKNjt5WKkp5eTk9evTg5ptvZtiwYRQWFnLs2LGYtr377rtZsGDBiddz585l4cKFDpX0K5EhCLEsiaKBRqkogUCALVu2cOutt1JaWkpGRgbz58+Padvp06fz3HPPAdYNiqKiIqZNm+ZkcQGdBUEp18nJySEnJ4dRo6ynV06aNIktW7bEtG2fPn04++yzKS0t5c9//jPDhg3j7LPPdrK4JyR7oNE2GqWinHvuufTq1Ytt27Zx0UUXsXbtWgYNGhTz9oWFhTz77LN8/vnnTJ8+3cGSfiWesyA4RQONUo08/vjjTJ06lfr6evr168eSJUti3vbaa6/l7rvvpqGhgRdeaHG6o7jQfjRKudDQoUNj7nPTmN/v54orrqBbt254ve1Xy0j2u04aaJSKo1AoxMaNG3n55ZfbbZ82p1tJCG0MVipOysrK6N+/P2PGjCE3N7fd9htpo4llSZTUq9F40uzP0zTXwTmO7M65hP2xSz4vSIFz+bdq3ii7c1nZPQ+tmY/LYYMGDeLTT5udDMAR2kaTCCkwr5PReZ1apvM6nUTbaJRSjnLDEAQNNEq5nPajUUo5zrrrpNOtKKUcpJdOSql2oYFGKeUobaNRSjlO+9EopRznhiEIGmiUcjm9dGonoVCIQCBAMBhMdFGUSohkv3QSY5rtxm2zH3n7OHLkCO+++y6ZmZlUV1fj8Xjwer34fD7GjptgDUOIlccHoYBz6b0+CNpID9bgpUDsQVN8Hkwg5Fj+to/B7mfUmm1akX7EsK/HPK9TMBgkPz+f7OxsXnvttdj30zptGh/RNT/XXFKyIKa0a+XqzcaY/LbsrzWSOww24ciRI4RCIbKysjh06NDJK0MN8BMb8fFx+/MD8YiN9HfanAcKrAGJNsYWmRE2547KTrc/dsnuXFZ2PiOwPie758HmebbzjJmFCxcycOBADh8+HPs+EsQN/Whc+ZiI7t2706FDB7p3757ooqgUVFVVxeuvv05hYWGiixIzfUyEUi5z++238+CDD3LkyJFEFyUmITxJPwTBlTUapZzy2muv0bNnT0aMGJHootgSz1kQRGS8iGwTkR0iMqeZdBNFxIhIi20+WqNRKsqGDRtYuXIlq1atoq6ujsOHDzNt2jSef/75RBetSfFsoxERL/AEMBaoAt4XkZXGmLJG6boAs4D3YslXazRKRXnggQeoqqqioqKCoqIiRo8endRBBqxbw3Fso7kE2GGM+dQYUw8UAdecJt2vgf8AYroLoYFGKdezhiDEsgCZIlIStTR+pms2UBn1uir83ld7ExkO9DLGvB5rCfXSSakmXH755Vx++eWJLkaLbF46fdmWfjQi4gF+B9xkZzsNNEq5nEE4Hr+xTtVAr6jXOeH3IroAecA6EQE4F1gpIgXGmCY7KmmgUcrl4jx6+30gV0T6YgWYycCPTuzLmENAZuS1iKwDZjcXZEADjVIpIV53nYwxARG5DXgT8ALPGGM+EpH7gBJjzMrW5KuBRimXi/cQBGPMKmBVo/fubiLt5bHkmXqBxpNmjV+KOX0rJiK708EJ5yLbjHBwkjq7+ds9BrufUWQbO+dBfPbPc4oyCMFQco91Sr1PP9QAhTYG2z1lf3Ceo4MwwfqS2p2kzu6gR7v52x1IauczglYNkrR9nlOUCQnH65J7CELqBRqlzjDGCMGA1miUUk4yaKBRSjnLGCHQoIFGKeUoIRRM7q9ycpdOKdUyA+ilk1LKUSGBuuT+Kid36ZRSsbH5LPj2llKBRqdbUWck64E0SS0lnkdjjKGuro66OhszAagzTn5+PuPHj090MeIvEmhiWRLE1TUaYwwNDQ0cP36cDh060KFDcveOVIkVy3QrlZWV/PjHP2bPnj2ICLfccguzZs1qh9K1gQFsTGWWCK6cQK6uro4NGzbQoUMHDh8+THp6OuFnYzD6ynH2JpATHxhnJy5zfDI1pyd4a49jtnseWpHexPB3sXv3bnbv3s3w4cM5cuQII0aM4JVXXmHQoEGx78u+No2PkAH5hsUxzln1L6ITyMVq79691NbWkpeXR1lZ2ckrQw1wg434uEzgZhvplzg8NgqscT/zbGxzj8BvbaSf24r8nZzcDazPye55sHueY5CVlUVWVhYAXbp0YeDAgVRXVzsdaNrGBW00rgw03bt3JyMjg27duiW6KCqFVVRUUFpayqhRoxJdlOZpoHGGz+fKYisXOXr0KBMnTmTBggV07do10cVpngYapdynoaGBiRMnMnXqVK677rpEF6dlGmiUchdjDDNmzGDgwIHccccdiS5O7JI80KREPxql4mXDhg0sW7aMt956i6FDhzJ06FBWrVrV8oaJFMKaxi2WJUG0RqNUlG9961u00OUj+eilk1LKcRpolFKO00CjlGoXGmiUUo7SGk0CeNJi7m4OWGNmljg8n5Cd+Yoi29xjcx9zbaa3m7+Tcy5FtrF7Huye51QVAmoTXYjmpd6nH2qAKTbuGixvxZgZJ+eNAutLOtvGNg+3YiyS3fydnHMJrM/J7nmwe55TlQGS/FFMqRdolDoT6aWTUspR2kajlHKcBhqllOMiQxCSmAYapVKB1miUUo7SSyellONc8HDylHlMhDGGmpqaRBdDJbGUnm4lGOOSICkRaCJBRh/xqZpTUlLC6tWrW0y3evVqLrroIvr378/8+fPboWRtFOd5nURkvIhsE5EdIjLnNOvvEJEyEfm7iKwVkfNbytP1gSYSZPx+P2lpaYkujnK5YDDIzJkzeeONNygrK2P58uWnzrSRbAzWEIRYlhaIiBd4AvgOMAiYIiKNp4AoBfKNMRcDK4AHW8rX1VWA+vp6ampq6NChAz6fD2MMnbt046id7uZ2x0Z50qzu8nbS2x3340mzuv3bSW9nLFJr8rc1vsvmZxTZxu55sHmejxw5QpcuXZpNtmnTJvr370+/fv0AmDx5MsXFxck/3Ur8LosuAXYYYz4FEJEi4BrgRLQ1xvxPVPqNwLSWMnVtoAmFQpSUlJwUZLxeLyteLqKhoYGOHTuemFQunhoaGggGg6Snp8c9b7CCp4g4Vjurq6sjLS0Nr9frSP41NTWkp6fj8cS/shypvXbs2NFW/nPmzGHfvn3079+f7Oxsevbs2eQlVHV1Nb169TrxOicnh/fee6/NZXeUvbtOmSISPdvcYmPM4qjX2UBl1OsqoLn5ZmYAb7S0U1deOtXW1lJTU8OAAQNOCjLHjx93NMgYY6ivr3d06t1AIOBYEAAQEUcfVen3+6mvr3ckbxEhPT2d2tpaW8cwf/58/vCHP5CVlUV1dTV79+5NrUZhe200Xxpj8qOWxafNMwYiMg3IBx5qKa0rA82RI0fo2LEj3bt3PynIHD9+3LEgA3D8+HH8fr9j+RtjMMY4UhuIEBFCoZBj+ft8PoLBoGPBzOv14vf7qauz3xX2d7/7XYvBJjs7m8rKr37Qq6qqyM7OblOZHRe5vR3L0rJqoFfU65zweycRkSuBuUCBMeZ4S5m6MtD07NkTr9dLMBg8Kch06tTJsSAQDAYJhUKONjhHjsdJHo/H8Ydvp6Wl0dDgXMeOtLQ0RKRVNaeWgs3IkSPZvn075eXl1NfXU1RUREFBQbyK7pz43d5+H8gVkb4i4gcmAyujE4jIMOD3WEFmbyyZurKNJvLLHwgE6NevHx9//DGXXnqpY5c0oVCI0tJShgwZQqdOnRzZB0B5eTkZGRn07NnTsX0cPHiQvXv3cuGFFzq2j0AgQGlpKSNGjHA08H/wwQdccMEFnHXWWba2ffvtt1m+fDkvvvgiL7/88knrfD4fixYtYty4cQSDQaZPn87gwYPjWfT4i+NYJ2NMQERuA94EvMAzxpiPROQ+oMQYsxLrUqkz8HL4/O40xjQbjaWFX7eknHeirq6OjRs3EgqFCAQCeDweRy83IkHN6dvnkfYZp76cYB1LMBh0vM9RexxLKBQiFAq16ljmzJnDF198gTGGzMxMMjMzY+pj45A2fUiSkW8YVNJyQoAS2WyMyW/L/lrDlYFGqRTTtkDTKd8wIMZAU5qYQOPKSyelVCM6qFIp5Sgdva2Ucpw++Eop5Tit0Sil2oUGGqWUo1zw4CsNNEq5nQsmkHPlEASl4mn//v2MHTuW3Nxcxo4dy4EDB06bbunSpeTm5pKbm8vSpUtPWV9QUEBeXp7TxT1VnB985QQNNOqMN3/+fMaMGcP27dsZM2bMaZ+qt3//fubNm8d7773Hpk2bmDdv3kkB6Y9//COdO3duz2J/JTL3dhwefOUUDTTqjFdcXExBQQFjx45l8eLFLFq06JRazZtvvsnYsWN59dVXueSSS9i3bx+//OUvAdi7dy+FhYVs3LiRHTt2MGfOKU+/dJ4+M1ip5LZnzx7uuece/va3v+HxeDh+/PgptZrPPvuM9evXU1hYSNeuXZk6dSovvfQSBw4c4Jprrjkxpqtnz55s2LCBN95o8VlQ8WViXBJEA406I1x55ZXk5eWdshQXF2OMobi4mFdffZWysjKMMbz44osnbb9p0yaOHz/O9OnT+dnPfsb69evJzc1l3rx5lJWVsWvXLpYuXcrnn3/O0KFDqaqqStCRJicdVKnOeL1792bXrl289tprzJw5k4qKCkSEQOCr1tOLL76YhoYGfD4ftbW1fPrpp4wbN46dO3dSXl5+4nlFgUAAn8/Htm3bTjx3OAZtG1Qp+QZiHFRJYgZVao1GnTGaqtWce+65GGOYOXMmEydOZNiwYQSDQXr37n3iEqqhoYGKigrS0tL461//ijGGtWvXsnPnTm699VbKy8v5xz/+gdfrxev1cvHFF/Pwww8n+IiTh/ajUWeMNWvWnPb9JUuWUFhYyK5du9i8eTPl5eX4/X4mTJjAQw89REFBAV6vl27durFt2zaGDBlCeno6xhj69OnDhg0beOSRR7j++uvxeDwEAgG++93vtuORRW47JS+t0agz3oABAzj77LPJy8vjN7/5DSJCnz598Pl8zJ49m+LiYrKzs+nQoQPdu3cnLy+PtLQ0gsEgV155JR9++CFz584lGAxy3nnnkZmZyZAhQ9rxCOL70GAnaKBRZ7yRI0dijGHnzp1ce+21HDp0iB/+8IeAdcm0bNkyCgoKOHz4MLNmzWL9+vUEAgG+9rWvMXHiRBoaGrj//vvZuHEjn332GV27dmXLli3teATJ32NPA4064/l8PubMmcPRo0cxxvD1r3+djh07UlZWxoEDBxg9ejQzZszA7/ezYMEC/H4/q1evJhQK8c1vfpPOnTuTnZ3N559/zrRp0/j4448ZPnx4Ox5B8tdotI1GKWDWrFk88cQTPPTQQzz55JMUFRXxwgsvsHKlNQFAeno6v/rVr/jwww8pLy9n3bp1jB49mmAwiIjQvXt3unbtyvr16+nTpw8HDx7E4/GQnp7Obbfd5nDpk39UpQYapfhq9oNZs2ZRUVHBT3/6U3Jzc1mwYAG/+tWvAJgxYwY33HADW7ZsobS0lI0bN7JixQqysrL45JNPyM3Nxe/3A/CDH/yAc845px2CDHw1+Xby0n40SjWyatUqbr/99hPTrcydO5e7776b/Px8CgoKqKur44YbbqC0tJTu3btTVFR0Sp+Ze++9l86dOzN79uxYdtnGfjR5Bl5uOSEAg3QWBKXOUG0MNIMNLI8x9dd1FgSlVGsk/7M8NdAo5XrJ3xist7eVagdPPvkkQ4cOZejQofTt25crrrgijrknfz8abaNRqh01NDQwevRofv7zn/O9730v8nYb22guNPCfMaYeq200SqW6WbNmMXr06OggEwfJf+mkgUapdvLss8/y2WefsWjRIgdy18Zgpc54mzdv5uGHH2b9+vV4PPFuGtUajVIKWLRoEfv37z/RCJyfn89TTz0Vp9w10CilsJ554xztR6OUclzyP/hKA41SrqeXTkopxyX/pZP2DFbK9eL74CsRGS8i20Rkh4icMhueiHQQkRfD698TkT4t5amBRinXi98QBBHxAk8A3wEGAVNEZFCjZDOAA8aY/sCjwH+0lK8GGqVcL66Tb18C7DDGfGqMqQeKgGsapbkGWBr+/wpgjESm6mxCS200bRqDoZRqD7vfhHszY0ycLiLRs80tNsYsjnqdDVRGva4CRjXK40QaY0xARA4BZwNfNrVTbQxWyuWMMeMTXYaW6KWTUipaNdAr6nVO+L3TphERH3AWsK+5TDXQKKWivQ/kikhfEfEDk4GVjdKsBG4M/38S8JZp4XkzeumklDoh3OZyG/Am4AWeMcZ8JCL3ASXGmJXA08AyEdkB7McKRs1q6cFXSinVZnrppJRynAYapZTjNNAopRyngUYp5TgNNEopx2mgUUo5TgONUspx/x/F99/1BsmvLwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solving time step:  4\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARoAAAD9CAYAAABjhbg0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAp3klEQVR4nO3df3wU1b3w8c83+4NA+FUMSExAQKICQfkRBK99tQqiSNtQhVooWBVyvY+P3uLjQ/vQ2lqxvcrVWkHxXktVQH1JVG5rUBFb8HqlaMRIatW0CBggCb/kp0B+7u55/phdXAJJdpKd7M7yfb9e+yK7c+bMmRn2u2fOnDNHjDEopZST0hJdAKVU6tNAo5RynAYapZTjNNAopRyngUYp5TgNNEopx2mgaYGIHBeRQc0su0VE/tLRZVLKjRIeaETkByJSGv5S7xGRN0Tk6+3Iz4jI4CafdROR34rIDhE5ISK7RGSViIxtKS9jTFdjzOdtLIdfRO4Tka3hbe4QkWdEZEBb8nNKuIzPJ7ocThKRAeH/F95El+VsldBAIyJ3A4uAB4Bzgf7AfwBT2pDXGf8TiUgn4C1gOPBtoDswBCgCrrOTl02rgALgB0AP4FLgQ2BCHPJWcaZByGHGmIS8sL58x4HvtZDmMuA94AiwB1gC+KOWG+AOYCtQAbwT/uxEOO/vA4XhdTNaKc8peUV9Njj89znAauBLYBPwK+AvzeR1NVAL9Gthe+eF8zsEbAP+OWrZfcDLwPPAMeBj4ELgp8B+oBK4Jir928CD4XJ9CRQDvcLLrgSqmmx7R7iMk4AGoDF8vD6KOjdPh49bNfBrwNPMfnQGVgCHgb8DP4neXng//wv4InyOfhS1rBPWD83u8GsR0Cm63OH89ofL8l1gMvBZ+Lj9LCqvNGA+sB04CLwUdQx2hc/l8fDrcuAWYCPwaDj9A+E8h0fl2QeoAXon6nuSKq/Ebdj6Tx4AvC2kGQ2MA7zAgPB/5Luilhvgz0AvoHPUZ4Oj0hQBy2MoT4t5hfN5CcgA8sJfwOYCzULgf1rZ3jtYtbd0YET4izg+vOw+oA64Nrzvz4a/pPcAPuCfCQfDcPq3w+XJC5fvv4Dnw8uupJlAE7Wt55ss/yPwu3BefbAC2L+0tK/A14Ac4G+R7YW//B8C9wJ+YBDwOXBtePn9QEl4G72Bd4FfRZU7EF43ss9fAC8A3YBhWMF8YDj93HBeOVgB7HfAyvCyAeFz6Y0q9y3h/P81fIw7h8/Hv0elmQu8mugvaSq8ErdhmAnstbnOXcAfo96byJezyWfRgWYdsDDq/QisGtKXwJZY8gI8WL/6F0cte4DmA83vgaIW9qMfEAS6RX32IOGAGP7y/zlq2Xewfok94ffdwmXrGX7/dpN9HIpVU/FgM9BgXcLWEw624c9mAP/dzL6cDBzh94V8FWjGAruapP8psCz893ZgctSya4Ed4b+vxAokTfd5bFT6D4Hvhv/+OzAhallW+JxFfqTOFGialm0sVu1Hwu9LgRsT9R1JpVcir0sPApki4jXGBM6UQEQuBH4L5ANdsP7TfNgkWWUM28mKvDHG/BXoKSJXA0/FmFfv8Lajl+9sZZsXtrD8POCQMeZYk/zyo97vi/q7FjhgjAlGvQfoihU0OUPZfEBmC2VozvnhdfeISOSzNJo/Nuc1WRb99/nAeSJyJOozD7Ahat3o47gz/FnEwTPsc9Pj0jVqW38UkVDU8iBW4GzOKftkjHlfRGqAK0VkD9aPzOoW1lcxSmRj8HtYv5zfbSHNfwL/AHKNMd2BnwHSJE1rw8/XA9eISEYMZWoury+wqtn9oj7r30I+64DLRCSnmeW7gV4i0q1JftUxlLE5TcvWCBzAaq/qElkgIh6swBnRdJ8rsc5LpjGmZ/jV3RgzrJnt7sG6XDlTOSqxLvF6Rr26GWMmh5fvxgoQ0eXe3eJeNq8SuK7JttKNMdVn2MeIM32+ApgF3ASsMsbUtbE8KkrCAo0x5ijW9fcTIvJdEekiIj4RuU5EHgon64Z1iXNcRC4Gbo8h631YbQERz2J9Gf4oInki4hGRdE6tPbRW1iDwB+C+cDmHAje3kH4dVnvPH0VktIh4w7fY/5eIzDbGVGK1RzwoIukicgkwB6vxt61michQEemC1faxKlzuz4B0EfmWiPiAn2O1YUTsAwaISFq47HuAPwGPiEh3EUkTkQtE5JvNbPcl4Kci8jURyQbujFq2CTgmIv9PRDqHj32eiIwJL18J/FxEeotIJtb/h7YegyeBfxOR8wHCeUbuXn4BhDj1/0Vzngeuxwo2z7axLKqJhN7eNsY8AtyN9Z//C6xfpTuBV8JJ5mHdHj6G1e7xYgzZ3gesEJEjInJj+BfpKqAceJ1w2wwwBrjRRnHvxKqm7wWWA8taST8NWBMu81HgE6zgti68fAZW28FurMbXX4YDVFs9Fy7XXqwG5h/ByYD+v7EuE6uxajhVUeu9HP73oIhsDv/9Q6zG23Ksu0mriLr8bOL+cH4V4X1bhVUjigTob2O1i1Vg1bCewrqrBdbdrFKsBuSPgc3hz9piMdZlzp9E5BhWw/DYcDlqgH8DNob/X4xrLpPwj8BmrNrOhubSKXsijV7KxUTkbawG3aZtTokoy+3AdGNMczWgpCcizwC7jTE/T3RZUoV2UlLtIiJZWJck7wG5wP/F6u/kSuGe2zcAIxNclJSS8CEIyvX8WH1WjmH1wC7G6o/iOiLyK6xL3IeNMRWJLk8ihIfJ7BeRT5pZLiLymIhsE5G/iciomPLVSyelVISIfAOrz9azxpi8MyyfjNXJcTJWG9hiY0yLYwZBazRKqSjGmHewhmI0ZwpWEDLGmBKsPmnN3Sg4qbU2Gq3uKOW8pn3DbBksYmpiTLsHPsUa3hKx1Biz1Mbmsjm1o2NV+LM9La2kjcFKuVwN8C8xpr0P6owxMfchixcNNEq5nNChX+RqTu39nUMMPdq1jUYpl0vDGnoeyysOVgM/DN99GgccDfcmb5HWaJRyOcEaBRuXvERWYo2czxSRKuCXkeyNMU9i9XafjPUMpRrg1ljy1UCjlMvF89LJGDOjleWRB8TZooFGKZeLZ43GKRpolHK5Dm4MbpNkL59SqhVuqNHoXacU88EHH3DJJZdQV1fHiRMnGDZsGJ98csZhKypFdPBdpzbRGk2KGTNmDAUFBfz85z+ntraWWbNmkZd32pAVlULcUKNpbVClDkFwoYaGBsaMGUN6ejrvvvsuHo8n0UVSLWvXEIRcEfNYjGknw4faM1jFxcGDBzl+/DiNjY3U1dWRkRHL45KVW2mNRiVEQUEB06dPp6Kigj179rBkiWufQ3W2aFeN5iIR87sY016lNRoVD88++yw+n48f/OAHBINB/umf/om33nqL8ePHJ7poyiGRxuBkpjUapRKvXTWaoSIm1qkjRmuNRinVFtphTynlODc0BmugUcrltEajlHKc1miUUo4Tkv+ukwYapVxOAF+s3+SAkyVpngYapVxOBLwaaJRSThIBX5IPZ9NAo5TL2arRJEiSF08p1RoR8HVKdClapoFGKbdzQUealH3C3tq1a7nooosYPHgwCxcujFu+lZWVXHXVVQwdOpRhw4axePHiuOUdEQwGGTlyJN/+9rfjmu+RI0eYNm0aF198MUOGDOG9996LW96PPvoow4YNIy8vjxkzZlBXV9f6Si2YPXs2ffr0OeWhXYcOHWLixInk5uYyceJEDh8+3N5ip4ZIoInllSApGWiCwSB33HEHb7zxBuXl5axcuZLy8vK45O31ennkkUcoLy+npKSEJ554Im55RyxevJghQ4bENU+AuXPnMmnSJP7xj3/w0UcfxW0b1dXVPPbYY5SWlvLJJ58QDAYpKipqV5633HILa9euPeWzhQsXMmHCBLZu3cqECRPi+gPiehpoOt6mTZsYPHgwgwYNwu/3M336dIqLi+OSd1ZWFqNGjQKgW7duDBkyhOrqVmcEjVlVVRWvv/46hYWFccsT4OjRo7zzzjvMmTMHAL/fT8+ePeOWfyAQoLa2lkAgQE1NDeedd1678vvGN75Br169TvmsuLiYm2++GYCbb76ZV155pV3bSBkCeGJ8JUhKBprq6mr69ftqeuCcnJy4BoOIHTt2UFZWxtixY+OW51133cVDDz1EWlp8T01FRQW9e/fm1ltvZeTIkRQWFnLixIm45J2dnc28efPo378/WVlZ9OjRg2uuuSYueUfbt28fWVlZAPTt25d9+/bFfRuupJdOqev48eNMnTqVRYsW0b1797jk+dprr9GnTx9Gjx4dl/yiBQIBNm/ezO23305ZWRkZGRlxu/Q4fPgwxcXFVFRUsHv3bk6cOMHzz8f6hJS2ERFE2vUYl9QhQKcYXwniykATCoVo6YFd2dnZVFZWnnxfVVVFdnZ23Lbf2NjI1KlTmTlzJjfccEPc8t24cSOrV69mwIABTJ8+nbfeeotZs2bFJe+cnBxycnJO1r6mTZvG5s2b45L3unXrGDhwIL1798bn83HDDTfw7rvvxiXvaOeeey579ljzye/Zs4c+ffrEfRuupDUaZzQ0NNDQ0NBssBkzZgxbt26loqKChoYGioqKKCgoiMu2jTHMmTOHIUOGcPfdd8clz4gHH3yQqqoqduzYQVFREePHj49bzaBv377069ePLVu2ALB+/XqGDh0al7z79+9PSUkJNTU1GGNYv369I43ZBQUFrFixAoAVK1YwZcqUuG/DlVwQaJL87nvzdu3aRf/+/fH7/adVob1eL0uWLOHaa68lGAwye/Zshg0bFpftbty4keeee47hw4czYsQIAB544AEmT54cl/yd9PjjjzNz5kwaGhoYNGgQy5Yti0u+Y8eOZdq0aYwaNQqv18vIkSO57bbb2pXnjBkzePvttzlw4AA5OTksWLCA+fPnc+ONN/L0009z/vnn89JLL8Wl/CkhyYcguPKZwXV1dbz33nuMGzeO6upqLrjgAr1eV27Wrv+8+d3ElMbYrCf/o88Mti0tLe1kW4wGG3XW0p7BHaOyspLt27ef1mYzadIkx7bp1rydzt+tebuaC+46JXkcjN2ZajYHDhxwbHtuzdvp/N2at6u5oEaT5MWzp7KykouH5RFsqD/5WauXUx4vBG08DSgqfUyXanbzD69j6zKwDfvgZP7i9bVpn2NZ52S57e6z14dpbLBXJrfQQNPxgg31UGyjDXuKwOs20n9LYL2N9BMESmy2qY8T+MjGOpd2QHo7+zDO5jEC6zjZPQ92z3OqigxBSGIp0Uaj1Fktzv1oRGSSiGwRkW0iMv8My/uLyH+LSJmI/E1EWu3bkXI1GqXOOpHG4HhkJeIBngAmAlXAByKy2hgT/YiCnwMvGWP+U0SGAmuAAS3lqzUapdwuvjWay4BtxpjPjTENQBHQtAu2ASID/HoAu1vLVGs0SrmdvcbgTBEpjXq/1BizNOp9NlAZ9b4KaPp4gvuAP4nIvwIZwNWtbVQDjVKpIPZv8oE49AyeASw3xjwiIpcDz4lInjEm1P7iKaWSU3zvOlUD/aLe54Q/izYHmARgjHlPRNKBTGB/c5lqG41SbhffNpoPgFwRGSgifmA6sLpJml3ABAARGQKkA1+0lKkrazQHDx4kEEjQlHtKJZs43nUyxgRE5E7gTax60jPGmE9F5H6g1BizGvi/wO9F5P9gNQzfYloZne3KQNOlSxfq6+vZuXNnoouiVOLFuWewMWYN1i3r6M/ujfq7HLjCTp6uDDSdO3emS5cuHDt2jNraWtLT03Xktjp76RAE54gIeXl5HDhwgJqaGjp37mw90Nvjs9fd3OO1urPbST/BZvpxNoOgx2t1+0+m9Hb2we4xiqxj9zzYPc+pSgON8/x+Px6Ph5qaGjp16gTBRufHLm20kf4Km+OKwAoC25u9U3i6C9JgZ2Ps6c/32c/f7tgoO8cIrONk9zzYPc+pLMnHOrk+0AB4PB4yMjKora1NdFGU6nhao+k4IkKXLl0SXQyVxPLz88nMzDxtBkzXi+NdJ6doPxp11igtLW02yLh6rm8XzIKggUYpXD7XtwYapdzB1XN9u2Du7ZRpo1Eq3lwz17c2BiuVGpJ6rm/BGm2UxPTSSalmuGaubxdcOmmgUaoZrpnrWxuDlXKHGTNmcPnll7NlyxZycnJ4+umnmT9/Pn/+85/Jzc1l3bp1zJ9/2nO6k0eSB5rUa6Px+Jwfu3SFg+OKIutcYOM3wOu1hhXYSW8n/7aMjbJzjCLr2D0Pds9zC1auXHnGz9evXx/7NhLFBdOtpF6gCTY6P3bpUxvphwlpe4/Hnh4I9e1Kj/o9Mac/2imLLPN5zOn3yCDb+dvZh1DfrvaOEcCwNpwHu+c5VeldJ6WU41wwBEEDjVJupzUapZTjNNAopRyngUYp1SH0rpNSylFao1FKOU7vOimlHKc1GmcEg8FEF0Gp5KGBxhmHDh3i+PHj7NixA2NM8g7fV6ojaKBxRu/evcnIyADgxIkTeL1e/H6/Na+TUmcho3ednCEiDBgwgN27d9PY2Ehtba1Vs/H4nB8kOcxGeq/XGvtjh9fL0U5ZttLvkUGO5m9rH+weo8g6ds+D3fOcokwaNCT5g69S4uj7fD58Pp/VdhNshBIbg+3G2Zzg7VJ7gyRDfbvaGvAI1qDHy8z/xJx+k3yT68x/xZz+DZlqO3+7gzbbMpDU9iR1Ns9zqk63YgQCnlhr8zYmDoyjlAg0ER5PktcfVUKVlpbGlO7RRx/lqaeeQkQYPnw4y5YtIz09easMRoSgN9avcoOjZWmONmooFaW6uprHHnuM0tJSPvnkE4LBIEVFRYkuVquCHk9Mr0RJqRqNUvEQCASora3F5/NRU1PDeeedl+gitcggBJN8DILWaJSKkp2dzbx58+jfvz9ZWVn06NGDa665JtHFapFBCOCJ6ZUoGmiUinL48GGKi4upqKhg9+7dnDhxgueffz7RxWqRQWigU0yvRNFAo1SUdevWMXDgQHr37o3P5+OGG27g3XffTXSxWhS5dIrllSgaaJSK0r9/f0pKSqipqcEYw/r16xkyZEiii9WqeAYaEZkkIltEZJuInHHqBxG5UUTKReRTEXmhtTy1MVipKGPHjmXatGmMGjUKr9fLyJEjue222xJdrBZF2mjiQUQ8wBPARKAK+EBEVhtjyqPS5AI/Ba4wxhwWkVZn1tNAo1QTCxYsYMGCBYkuRsysS6e4fZUvA7YZY/XQFJEiYApQHpXmn4EnjDGHAYwx+1vLVAONUi5nNQb7Y02eKSLRPReXGmOWRr3PBiqj3lcBY5vkcSGAiGzEerbffcaYFrtbp16g8fqsYQWxsjs5mt1xP3bHIQF4PWySb8acXLwe3pCpjuWP12N7LJXt8V1tmaTO7nlOUQbsXDodMMbkt3OTXiAXuBLIAd4RkeHGmCMtrZBaAo32x8xstzH+44I025Ov2RlXBNbYon81D8Wc/nH5CSvNd2NOP0NesZ2/3bFRdo4RWMfJ7nmwfZ5TVlwvnaqBflHvc8KfRasC3jfGNAIVIvIZVuD5oLlM9a6TUi4X59vbHwC5IjJQRPzAdGB1kzSvYNVmEJFMrEupFkfdpl6NRqmzULz6yBhjAiJyJ/AmVvvLM8aYT0XkfqDUGLM6vOwaESkHgsCPjTEHW8pXA41SLhfvsU7GmDXAmiaf3Rv1twHuDr9iooFGKZczCPVJPg2CBhqlXM4No7c10CjlchpolFIdIpGPgIiFBhqlXC7OQxAckdylU0q1Si+dHHL48GFqa2vZtm0bDQ0NiAhpaWk6kZw6K1l3nWIe65QQrgw0GRkZ+Hw+unXrBkAoFCIYDBIKhayxTnbHzFxgo4O07TmRbI4rAsSbxuPyk5jTe7zCDHkl5vRpNvO3PzbK5jEKr2PrPLRhbFTKTreil07O8Pv9eL1ezj33XLZv337qwraMddrZGHv683225ziyM+cSWPMu2R27ZG6JPX9ZHrKdv915o9oyl5Xd82D3PMc63cqRI0coLCzkk08+QUR45plnuPzyy2PfVgLopZNSLjN37lwmTZrEqlWraGhooKamJtFFapG20SjlMkePHuWdd95h+fLlgFV79vuTu/3DDYFGR28rFaWiooLevXtz6623MnLkSAoLCzlx4kRM6957770sWrTo5Pt77rmHxYsXO1TSr0SGIMTyShQNNEpFCQQCbN68mdtvv52ysjIyMjJYuHBhTOvOnj2bZ599FrBuUBQVFTFr1iwniwvoLAhKuU5OTg45OTmMHWs9vXLatGls3rw5pnUHDBjAOeecQ1lZGX/6058YOXIk55xzjpPFPSnZA4220SgVpW/fvvTr148tW7Zw0UUXsX79eoYOHRrz+oWFhSxfvpy9e/cye/ZsB0v6lXjOguAUDTRKNfH4448zc+ZMGhoaGDRoEMuWLYt53euvv557772XxsZGXnih1emO4kL70SjlQiNGjIi5z01Tfr+fq666ip49e+LxdFwtI9nvOmmgUSqOQqEQJSUlvPzyyx22TZvTrSSENgYrFSfl5eUMHjyYCRMmkJub22HbjbTRxPJKlNSr0dgd6+T1Wt3ZbaS3M8eR7TmXgDSbY5e8ArI89vztjo1qy7xR9ueysnke2jIPlMOGDh3K55/bG3oRD9pGkwg6r1OrdF6n1KNtNEopR7lhCIIGGqVcTvvRKKUcZ9110ulWlFIO0ksnpVSH0ECjlHKUttEopRyn/WiUUo5zwxAEDTRKuZxeOnWQUChEIBAgGAwmuihKJUSyXzqJMS1247bRx7vjHDt2jPfee4/MzEyqq6tJS0vD4/Hg9XqZeN1kaxhCrDxeCAZiT+/1QsBOeg8E7AVA8aZhArF3x0/zCqFA7KfKbv6298HuMWrLOnbPm8fL6BGXxjyvUzAYJD8/n+zsbF577bXYt9M27Rof0T0/11xWuiimtOvl2x8aY/Lbs722SO4w2Ixjx44RCoXIysri6NGjpy4MNEKJjfg4TmyPmUnbezzm5KG+Xds0x5HdsUV2512ym7/duazsHCOwjpPtsUs2z7OdZ8wsXryYIUOG8OWXX8a+jQRxQz8aVz4molevXnTq1IlevXoluigqBVVVVfH6669TWFiY6KLETB8ToZTL3HXXXTz00EMcO3Ys0UWJSYi0pB+C4MoajVJOee211+jTpw+jR49OdFFsiecsCCIySUS2iMg2EZnfQrqpImJEpNU2H63RKBVl48aNrF69mjVr1lBXV8eXX37JrFmzeP755xNdtGbFs41GRDzAE8BEoAr4QERWG2PKm6TrBswF3o8lX63RKBXlwQcfpKqqih07dlBUVMT48eOTOsiAdWs4jm00lwHbjDGfG2MagCJgyhnS/Qr4d6Aulkw10CjletYQhFheQKaIlEa9bmuSWTZQGfW+KvzZV1sTGQX0M8a8HmsJ9dJJqWZceeWVXHnllYkuRqtsXjodaE8/GhFJA34L3GJnPQ00SrmcQaiP31inaqBf1Puc8GcR3YA84G0RAegLrBaRAmNMsx2VNNAo5XJxHr39AZArIgOxAsx04Acnt2XMUSAz8l5E3gbmtRRkQAONUikhXnedjDEBEbkTeBPwAM8YYz4VkfuBUmPM6rbkq4FGKZeL9xAEY8waYE2Tz+5tJu2VseSZeoHG67PGL8XK7kRkXq81LsdGevuTqXnYJN+MOXlbJnizk7/tCeHsHiNo24Rwds9zijIIwVByj3VKvaMfaIT1NgbbTRDYaCP9FQKf2kg/zN4gTLAGGNqdpM7uoEe7+dsdSGrrGAEMa8N5sHueU5QJCfV1yT0EIfUCjVJnGWOEYEBrNEopJxk00CilnGWMEGjUQKOUcpQQCib3Vzm5S6eUap0B9NJJKeWokEBdcn+Vk7t0SqnY2HwWfEdLqUCj062os5L1QJqklhLPozHGUFdXR11dTM/gUWep/Px8Jk2alOhixF8k0MTyShBX12iMMTQ2NlJfX0+nTp3o1Cm5e0eqxIplupXKykp++MMfsm/fPkSE2267jblz53ZA6drBADamMksE1waaUChEaWkpgUCAjIwMws/GsMY62elu7vFa3dntpB/m4Nio8DpHO2XZSm93LJLd/G3tg91jFFnH7nmwe55j4PV6eeSRRxg1ahTHjh1j9OjRTJw4kaFDh8a+rY5mgPpEF6Jlrgw0+/fvp7a2lry8PMrLy09dGGiE122MgflWG8bM2B2TY2diNLAGF263MZPkBWmw08ZP2vk++/nbndzNzjGCto1dsnueY5CVlUVWlhWEu3XrxpAhQ6iurk7+QJPkbTSuDDS9evUiIyODnj17JrooKoXt2LGDsrIyxo4dm+iitEwDjTO8XlcWW7nI8ePHmTp1KosWLaJ79+6JLk7LNNAo5T6NjY1MnTqVmTNncsMNNyS6OK3TQKOUuxhjmDNnDkOGDOHuu+9OdHFil+SBJiX60SgVLxs3buS5557jrbfeYsSIEYwYMYI1a9a0vmIihbCmcYvllSBao1Eqyte//nWMsXnHLNH00kkp5TgNNEopx2mgUUp1CA00SilHaY0mATy+mLubW+nbMGbG7pgcO/MVRda5wMYNQa/XGlZgJ72d/Nsy55KdYxRZx+55sHueU1UIqE10IVqWekc/2AjFNu4aTGnDmBm7Y3JKbN7FGGdzfNSlHZDezj6Ms3mMoG1jl+ye51RlgCR/FFPqBRqlzkZ66aSUcpS20SilHKeBRinluMgQhCSmgUapVKA1GqWUo/TSSSnlOBc8nDxlHhNhjKGmpibRxVBJLKWnWwnG+EqQlAg0kSCjj/hULSktLWXt2rWtplu7di0XXXQRgwcPZuHChR1QsnaK87xOIjJJRLaIyDYRmX+G5XeLSLmI/E1E1ovI+a3l6fpAEwkyfr8fn89GN3ylziAYDHLHHXfwxhtvUF5ezsqVK0+faSPZGKwhCLG8WiEiHuAJ4DpgKDBDRJpOAVEG5BtjLgFWAQ+1lq+rqwANDQ3U1NTQqVMnvF4vxhi69ujJcTvdze2OjbI7b5TXZ3XJt8Prsze2qCPS29kHu8cI2jBGzWdvWIHHx7Fjx+jWrVuLyTZt2sTgwYMZNMiaJ2v69OkUFxcn/3Qr8bssugzYZoz5HEBEioApwMloa4z576j0JcCs1jJ1baCJTCAXHWQ8Hg+rXiyisbGRzp07fzWpXBw1NjYSDAZJT0+Pe95gBU8Rcax2VldXh8/nw+PxOJJ/TU0N6enppKXFv7Icqb127tzZVv7z58/n4MGDDB48mOzsbPr06dPsJVR1dTX9+vU7+T4nJ4f333+/3WV3lL27TpkiEj1l51JjzNKo99lAZdT7KqCl+WbmAG+0tlFXXjrV1tZSU1PDxRdffEqQqa+vdzTIGGNoaGhwdOrdQCDgWBAAEBFHH1Xp9/tpaGhwJG8RIT09ndraWlv7sHDhQn7/+9+TlZVFdXU1+/fvT61GYXttNAeMMflRr6VnzDMGIjILyAcebi2tKwPNsWPH6Ny5M7169TolyNTX1zsWZADq6+vx+/2O5W+MwRjjSG0gQkQIhWzMUmmT1+slGAw6Fsw8Hg9+v5+6OvtdYX/729+2Gmyys7OprPzqB72qqors7Ox2ldlxkdvbsbxaVw30i3qfE/7sFCJyNXAPUGCMaXVCXlcGmj59+uDxeAgGg6cEmS5dujgWBILBIKFQyNEG58j+OCktLc3xh2/7fD4aG53r2OHz+RCRNtWcWgs2Y8aMYevWrVRUVNDQ0EBRUREFBQXxKrpz4nd7+wMgV0QGiogfmA6sjk4gIiOB32EFmf2xZOrKNprIL38gEGDQoEF89tlnjBs3zrFLmlAoRFlZGcOHD6dLly6ObAOgoqKCjIwM+vTp49g2jhw5wv79+7nwwgsd20YgEKCsrIzRo0c7Gvg/+ugjLrjgAnr06GFr3XfeeYeVK1fy4osv8vLLL5+yzOv1smTJEq699lqCwSCzZ89m2LBh8Sx6/MVxrJMxJiAidwJvAh7gGWPMpyJyP1BqjFmNdanUFXg5fH53GWNajMbSyq9bUs47UVdXR0lJCaFQiEAgQFpamqOXG5Gg5vTt80j7jFNfTrD2JRgMOt7nqCP2JRQKEQqF2rQv8+fP54svvsAYQ2ZmJpmZmTH1sXFIuw6SZOQbhpa2nhCgVD40xuS3Z3tt4cpAo1SKaV+g6ZJvuDjGQFOWmEDjyksnpVQTOqhSKeUoHb2tlHKcPvhKKeU4rdEopTqEBhqllKNc8OArDTRKuZ0LJpBz5RAEpeLp0KFDTJw4kdzcXCZOnMjhw4fPmG7FihXk5uaSm5vLihUrTlteUFBAXl6e08U9XZwffOUEDTTqrLdw4UImTJjA1q1bmTBhwhmfqnfo0CEWLFjA+++/z6ZNm1iwYMEpAekPf/gDXbt27chifyUy93YcHnzlFA006qxXXFxMQUEBEydOZOnSpSxZsuS0Ws2bb77JxIkTefXVV7nssss4ePAgP/vZzwDYv38/hYWFlJSUsG3bNubPP+3pl87TZwYrldz27dvHL3/5S/7617+SlpZGfX39abWanTt3smHDBgoLC+nevTszZ87kpZde4vDhw0yZMuXkmK4+ffqwceNG3nij1WdBxZeJ8ZUgGmjUWeHqq68mLy/vtFdxcTHGGIqLi3n11VcpLy/HGMOLL754yvqbNm2ivr6e2bNn8+Mf/5gNGzaQm5vLggULKC8vZ/fu3axYsYK9e/cyYsQIqqqqErSnyUkHVaqzXv/+/dm9ezevvfYad9xxBzt27EBECAS+aj295JJLaGxsxOv1Ultby+eff861117Lrl27qKioOPm8okAggNfrZcuWLSefOxyD9g2qlHwDMQ6qJDGDKrVGo84azdVq+vbtizGGO+64g6lTpzJy5EiCwSD9+/c/eQnV2NjIjh078Pl8/OUvf8EYw/r169m1axe33347FRUV/P3vf8fj8eDxeLjkkkv4zW9+k+A9Th7aj0adNdatW3fGz5ctW0ZhYSG7d+/mww8/pKKiAr/fz+TJk3n44YcpKCjA4/HQs2dPtmzZwvDhw0lPT8cYw4ABA9i4cSOPPPIIN954I2lpaQQCAb71rW914J5FbjslL63RqLPexRdfzDnnnENeXh6//vWvEREGDBiA1+tl3rx5FBcXk52dTadOnejVqxd5eXn4fD6CwSBXX301H3/8Mffccw/BYJDzzjuPzMxMhg8f3oF7EN+HBjtBA406640ZMwZjDLt27eL666/n6NGjfP/73wesS6bnnnuOgoICvvzyS+bOncuGDRsIBAJ87WtfY+rUqTQ2NvLAAw9QUlLCzp076d69O5s3b+7APUj+HnsaaNRZz+v1Mn/+fI4fP44xhksvvZTOnTtTXl7O4cOHGT9+PHPmzMHv97No0SL8fj9r164lFApxxRVX0LVrV7Kzs9m7dy+zZs3is88+Y9SoUR24B8lfo9E2GqWAuXPn8sQTT/Dwww/z5JNPUlRUxAsvvMDq1dYEAOnp6fziF7/g448/pqKigrfffpvx48cTDAYREXr16kX37t3ZsGEDAwYM4MiRI6SlpZGens6dd97pcOmTf1SlBhql+Gr2g7lz57Jjxw5+9KMfkZuby6JFi/jFL34BwJw5c7jpppvYvHkzZWVllJSUsGrVKrKysti+fTu5ubn4/X4Avve973Huued2QJCBrybfTl7aj0apJtasWcNdd911crqVe+65h3vvvZf8/HwKCgqoq6vjpptuoqysjF69elFUVHRan5n77ruPrl27Mm/evFg22c5+NHkGXm49IQBDdRYEpc5S7Qw0wwysjDH1pToLglKqLZL/WZ4aaJRyveRvDNbb20p1gCeffJIRI0YwYsQIBg4cyFVXXRXH3JO/H4220SjVgRobGxk/fjw/+clP+M53vhP5uJ1tNBca+I8YU0/UNhqlUt3cuXMZP358dJCJg+S/dNJAo1QHWb58OTt37mTJkiUO5K6NwUqd9T788EN+85vfsGHDBtLS4t00qjUapRSwZMkSDh06dLIROD8/n6eeeipOuWugUUphPfPGOdqPRinluOR/8JUGGqVcTy+dlFKOS/5LJ+0ZrJTrxffBVyIySUS2iMg2ETltNjwR6SQiL4aXvy8iA1rLUwONUq4XvyEIIuIBngCuA4YCM0RkaJNkc4DDxpjBwKPAv7eWrwYapVwvrpNvXwZsM8Z8boxpAIqAKU3STAFWhP9eBUyQyFSdzWitjaZdYzCUUh1hz5twX2aMidNFJHq2uaXGmKVR77OByqj3VcDYJnmcTGOMCYjIUeAc4EBzG9XGYKVczhgzKdFlaI1eOimlolUD/aLe54Q/O2MaEfECPYCDLWWqgUYpFe0DIFdEBoqIH5gOrG6SZjVwc/jvacBbppXnzeilk1LqpHCby53Am4AHeMYY86mI3A+UGmNWA08Dz4nINuAQVjBqUWsPvlJKqXbTSyellOM00CilHKeBRinlOA00SinHaaBRSjlOA41SynEaaJRSjvv/B2V6OaUyBbMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solving time step:  5\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARoAAAD9CAYAAABjhbg0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAp8UlEQVR4nO3de3hU1b3w8e8vMxnCVYqAxgQEJGogKJcA0vapCqJobahCLQhWhRzP6yun+PrSvrQqFdujHK0VLZ5a6gXUR6PS1iACWvBYKYIYiVaMRdFwSUCQ+yXXmaz3j72DQyTJ3snszOzh93meeWBmr1l77dmZ36y99v7tJcYYlFLKSynxboBSKvlpoFFKeU4DjVLKcxpolFKe00CjlPKcBhqllOc00DRBRI6KSL9Glt0kIv9o6zYp5UdxDzQicr2IFNlf6l0iskJEvtuK+oyI9G/wWmcR+Z2IbBWRYyKyXUSWiMjIpuoyxnQyxnzRwnaEROQeEfnMXudWEXlKRPq0pD6v2G18Lt7t8JKI9LH/LoLxbsupKq6BRkTuAOYD9wFnAL2B/wbGt6Cuk/4RiUg74E1gEHA10AXIBgqAK93U5dISIA+4HjgNuBB4HxgTg7pVjGkQ8pgxJi4PrC/fUeBHTZQZAawDDgK7gAVAKGq5AW4DPgNKgbft147Zdf8YyLff27GZ9pxQV9Rr/e3/nw4sBQ4DG4BfA/9opK7LgEqgVxPrO8uubz+wBfi3qGX3AC8DzwFHgI+Ac4FfAHuAHcDlUeXfAu6323UYKAS62csuAcoarHur3cZxQA1Qa39eH0btmyftz60c+A0QaGQ72gOLgQPAJ8DPo9dnb+efga/sffTTqGXtsH5odtqP+UC76Hbb9e2x2/JD4CrgU/tz+2VUXSnAbOBzYB/wUtRnsN3el0ftxyjgJmAt8LBd/j67zkFRdfYEKoAe8fqeJMsjfiu2/sjDQLCJMsOAi4Ag0Mf+Q749arkB/gZ0A9pHvdY/qkwBsMhBe5qsy67nJaAjkGN/ARsLNPOAvzezvrexem9pwGD7izjaXnYPUAVcYW/7M/aX9E4gFfg37GBol3/Lbk+O3b4/A8/Zyy6hkUATta7nGiz/K/BHu66eWAHs35vaVuBbQCbwz/r12V/+94E5QAjoB3wBXGEvvxdYb6+jB/AO8Ouodoft99Zv81fA80BnYCBWMO9rl59p15WJFcD+CLxgL+tj78tgVLtvsuv/D/szbm/vj/+KKjMTeDXeX9JkeMRvxTAF+NLle24H/hr13NR/ORu8Fh1oVgHzop4PxuohHQY2O6kLCGD96p8ftew+Gg80fwIKmtiOXkAE6Bz12v3YAdH+8v8tatkPsH6JA/bzznbbutrP32qwjQOweioBXAYarEPYauxga782GfifRrbleOCwn+fzdaAZCWxvUP4XwNP2/z8HropadgWw1f7/JViBpOE2j4wq/z7wQ/v/nwBjopal2/us/kfqZIGmYdtGYvV+xH5eBFwXr+9IMj3ieVy6D+guIkFjTPhkBUTkXOB3QC7QAeuP5v0GxXY4WE96/RNjzAdAVxG5DHjCYV097HVHL9/WzDrPbWL5WcB+Y8yRBvXlRj3fHfX/SmCvMSYS9RygE1bQ5CRtSwW6N9GGxpxtv3eXiNS/lkLjn81ZDZZF//9s4CwRORj1WgBYE/Xe6M9xm/1avX0n2eaGn0unqHX9VUTqopZHsAJnY07YJmPMuyJSAVwiIruwfmSWNvF+5VA8B4PXYf1y/rCJMn8A/gVkGWO6AL8EpEGZ5tLPVwOXi0hHB21qrK6vsLrZvaJe691EPauAESKS2cjynUA3EencoL5yB21sTMO21QJ7scarOtQvEJEAVuCs13Cbd2Dtl+7GmK72o4sxZmAj692FdbhysnbswDrE6xr16GyMucpevhMrQES3e2eTW9m4HcCVDdaVZowpP8k21jvZ64uBqcANwBJjTFUL26OixC3QGGMOYR1/PyYiPxSRDiKSKiJXisgDdrHOWIc4R0XkfOBWB1XvxhoLqPcM1pfhryKSIyIBEUnjxN5Dc22NAH8B7rHbOQC4sYnyq7DGe/4qIsNEJGifYv9fIjLNGLMDazzifhFJE5ELgOlYg78tNVVEBohIB6yxjyV2uz8F0kTk+yKSCtyFNYZRbzfQR0RS7LbvAt4AHhKRLiKSIiLniMjFjaz3JeAXIvItEckAZkQt2wAcEZH/JyLt7c8+R0SG28tfAO4SkR4i0h3r76Gln8HjwH+KyNkAdp31Zy+/Auo48e+iMc8B12AFm2da2BbVQFxPbxtjHgLuwPrj/wrrV2kG8IpdZBbW6eEjWOMeLzqo9h5gsYgcFJHr7F+kS4ES4DXssRlgOHCdi+bOwOqmfwksAp5upvxEYLnd5kPAJqzgtspePhlr7GAn1uDrr+wA1VLP2u36EmuA+adwPKD/b6zDxHKsHk5Z1Ptetv/dJyIb7f//BGvwtgTrbNISog4/G7jXrq/U3rYlWD2i+gB9Nda4WClWD+sJrLNaYJ3NKsIaQP4I2Gi/1hKPYB3mvCEiR7AGhkfa7agA/hNYa/9dXNRYJfaPwEas3s6axsopd+oHvZSPichbWAO6Dcec4tGWW4FJxpjGekAJT0SeAnYaY+6Kd1uShV6kpFpFRNKxDknWAVnA/8W63smX7Cu3rwWGxLkpSSXuKQjK90JY16wcwboCuxDrehTfEZFfYx3iPmiMKY13e+LBTpPZIyKbGlkuIvKoiGwRkX+KyFBH9eqhk1Kqnoh8D+uarWeMMTknWX4V1kWOV2GNgT1ijGkyZxC0R6OUimKMeRsrFaMx47GCkDHGrMe6Jq2xEwXHNTdGo90dpbzX8NowV/qLmAqHZXfBx1jpLfUWGmMWulhdBide6Fhmv7arqTfpYLBSPlcB/LvDsvdAlTHG8TVksaKBRimfE9r0i1zOiVd/Z+LginYdo1HK51KwUs+dPGJgKfAT++zTRcAh+2ryJmmPRimfE6ws2JjUJfICVuZ8dxEpA35VX70x5nGsq92vwrqHUgVws5N6NdAo5XOxPHQyxkxuZnn9DeJc0UCjlM/FskfjFQ00SvlcGw8Gt0iit08p1Qw/9Gj0rFOSee+997jggguoqqri2LFjDBw4kE2bTpq2opJEG591ahHt0SSZ4cOHk5eXx1133UVlZSVTp04lJ+cbKSsqifihR9NcUqWmIPhQTU0Nw4cPJy0tjXfeeYdAIBDvJqmmtSoFIUvEPOqw7FXwvl4ZrGJi3759HD16lNraWqqqqujY0cntkpVfaY9GxUVeXh6TJk2itLSUXbt2sWCBb+9DdapoVY/mPBHzR4dlL9UejYqFZ555htTUVK6//noikQjf/va3efPNNxk9enS8m6Y8Uj8YnMi0R6NU/LWqRzNAxDidOmKY9miUUi2hF+wppTznh8FgDTRK+Zz2aJRSntMejVLKc0Lin3XSQKOUzwmQ6vSbHPayJY3TQKOUz4lAUAONUspLIpCa4OlsGmiU8jlXPZo4SfDmKaWaIwKp7eLdiqZpoFHK73xwIU3S3mFv5cqVnHfeefTv35958+bFrN4dO3Zw6aWXMmDAAAYOHMgjjzwSs7rrRSIRhgwZwtVXXx3Teg8ePMjEiRM5//zzyc7OZt26dTGr++GHH2bgwIHk5OQwefJkqqqqmn9TE6ZNm0bPnj1PuGnX/v37GTt2LFlZWYwdO5YDBw60ttnJoT7QOHnESVIGmkgkwm233caKFSsoKSnhhRdeoKSkJCZ1B4NBHnroIUpKSli/fj2PPfZYzOqu98gjj5CdnR3TOgFmzpzJuHHj+Ne//sWHH34Ys3WUl5fz6KOPUlRUxKZNm4hEIhQUFLSqzptuuomVK1ee8Nq8efMYM2YMn332GWPGjInpD4jvaaBpexs2bKB///7069ePUCjEpEmTKCwsjEnd6enpDB06FIDOnTuTnZ1NeXmzM4I6VlZWxmuvvUZ+fn7M6gQ4dOgQb7/9NtOnTwcgFArRtWvXmNUfDoeprKwkHA5TUVHBWWed1ar6vve979GtW7cTXissLOTGG28E4MYbb+SVV15p1TqShgABh484ScpAU15eTq9eX08PnJmZGdNgUG/r1q0UFxczcuTImNV5++2388ADD5CSEttdU1paSo8ePbj55psZMmQI+fn5HDt2LCZ1Z2RkMGvWLHr37k16ejqnnXYal19+eUzqjrZ7927S09MBOPPMM9m9e3fM1+FLeuiUvI4ePcqECROYP38+Xbp0iUmdy5Yto2fPngwbNiwm9UULh8Ns3LiRW2+9leLiYjp27BizQ48DBw5QWFhIaWkpO3fu5NixYzz3nNM7pLSMiCDSqtu4JA8B2jl8xIkvA01dXR1N3bArIyODHTt2HH9eVlZGRkZGzNZfW1vLhAkTmDJlCtdee23M6l27di1Lly6lT58+TJo0iTfffJOpU6fGpO7MzEwyMzOP974mTpzIxo0bY1L3qlWr6Nu3Lz169CA1NZVrr72Wd955JyZ1RzvjjDPYtcuaT37Xrl307Nkz5uvwJe3ReKOmpoaamppGg83w4cP57LPPKC0tpaamhoKCAvLy8mKybmMM06dPJzs7mzvuuCMmdda7//77KSsrY+vWrRQUFDB69OiY9QzOPPNMevXqxebNmwFYvXo1AwYMiEndvXv3Zv369VRUVGCMYfXq1Z4MZufl5bF48WIAFi9ezPjx42O+Dl/yQaBJ8LPvjdu+fTu9e/cmFAp9owsdDAZZsGABV1xxBZFIhGnTpjFw4MCYrHft2rU8++yzDBo0iMGDBwNw3333cdVVV8Wkfi/9/ve/Z8qUKdTU1NCvXz+efvrpmNQ7cuRIJk6cyNChQwkGgwwZMoRbbrmlVXVOnjyZt956i71795KZmcncuXOZPXs21113HU8++SRnn302L730UkzanxQSPAXBl/cMrqqqYt26dVx00UWUl5dzzjnn6PG68rNW/fHmdhZT5HBYT/6u9wx2LSUl5fhYjAYbdcrSK4Pbxo4dO/j888+/MWYzbtw4z9bp17q9rt+vdfuaD846JXgcdO5kPZu9e/d6tj6/1u11/X6t29d80KNJ8Oa5s2PHDrJzcghXVx9/rdnDqWAAwhHnK4kq7+hQzW39gAQDrg4DJRjAuNwGV4eZLrdBUoOut9npOo63OxiEsPO7OElqkLqaWndt8gsNNG0vXF3N2eYTx+W3STbnmE2Oy38uOWQb59effCJDGWzcJS9+IKMYYf7uuPwGuZhR5k3H5dfJaNf1u9mGD2SUq88IrM/J7X5wu5+TVn0KQgJLijEapU5pMb6ORkTGichmEdkiIrNPsry3iPyPiBSLyD9FpNlrO5KuR6PUKad+MDgWVYkEgMeAsUAZ8J6ILDXGRN+i4C7gJWPMH0RkALAc6NNUvdqjUcrvYtujGQFsMcZ8YYypAQqAhpdgG6A+we80YGdzlWqPRim/czcY3F1EiqKeLzTGLIx6ngHsiHpeBjS8PcE9wBsi8h9AR+Cy5laqgUapZOD8m7w3BlcGTwYWGWMeEpFRwLMikmOMqWt985RSiSm2Z53KgV5RzzPt16JNB8YBGGPWiUga0B3Y01ilOkajlN/FdozmPSBLRPqKSAiYBCxtUGY7MAZARLKBNOCrpir1ZY9m3759hF1crKVUUovhWSdjTFhEZgCvY/WTnjLGfCwi9wJFxpilwP8F/iQi/wdrYPgm00x2ti8DTYcOHaiurmbbtm3xbopS8RfjK4ONMcuxTllHvzYn6v8lwHfc1OnLQNO+fXs6dOjAkSNHqKysJC0tTTO31alLUxC8IyLk5OSwd+9eKioqaN++PSkpKUhq0N3l5sEAn0tO8+Wiyn8iQ12V/0BGOS9vv2eDXOyq/DoZ7Wn9rrbB7Wdkv8ftfnC3n337p948DTTeC4VCBAIBKioqaNeuHaY27Dpn5lzzoePyn8qFrvN+3OQhgZWLNMYsc1x+tVzN5cb5dDJvyHjX9bvNpWpJfpfb/eB2Pye1BM918n2gAQgEAnTs2JHKysp4N0Wptqc9mrYjInTo0CHezVAJLDc3l+7du39jBkzfi+FZJ6/odTTqlFFUVNRokPH1XN8+mAVBA41S+Hyubw00SvmDr+f69sHc20kzRqNUrPlmrm8dDFYqOST0XN+ClW2UwPTQSalG+Gaubx8cOmmgUaoRvpnrWweDlfKHyZMnM2rUKDZv3kxmZiZPPvkks2fP5m9/+xtZWVmsWrWK2bO/cZ/uxJHggcbXc29/+9vf5p133jlh2ZgrLsfUuriFRCvmdfKkPO7naUq08i3Z5rbYD67+LtpW6+bePldM0QKHK7pC596OCVMb9nTeJbe5S+tkNFeaPzsuD7BCJnC9edJx+edlOjeZPzguv0hudV2/m21YIRNalN/l5dxRrpM8/UTPOimlPOeDFAQNNEr5nfZolFKe00CjlPKcBhqlVJvQG18ppTylPRqllOf0rJNSynPao/FGJOLyqlOlkpkGGm/s37+fo0ePsnXrVowxiZu+r1Rb0EDjjR49etCxY0cAjh07RjAYJBQKkZKiOaLq1GQS/KxTUiRV1tbWUlNTg4jwg2t+mFBJla4TEgEJpmDCdQlU3uMkTNCkylYYNlTMO/9wVjatoyZVtlhqaiqpqalEIhFMbTihkiRXyASmG4eptbYnZQZ3m186Lv9ruY/55hbH5W+Xha7rd7MNT8qMFiWSejlJ3QcyKmmnWzEC4YDT3rzzH5hYSopAUy8QSPD+o4qroqIiR+UefvhhnnjiCUSEQYMG8fTTT5OWlrj3yjQiRBxP+VvjaVsao4MaSkUpLy/n0UcfpaioiE2bNhGJRCgoKIh3s5oVCQQcPeIlqXo0SsVCOBymsrKS1NRUKioqOOuss+LdpCYZhEiC5yBoj0apKBkZGcyaNYvevXuTnp7OaaedxuWXXx7vZjXJIIQJOHrEiwYapaIcOHCAwsJCSktL2blzJ8eOHeO5556Ld7OaZBBqaOfoES8aaJSKsmrVKvr27UuPHj1ITU3l2muv/cZ9qRNN/aGTk0e8aKBRKkrv3r1Zv349FRUVGGNYvXo12dnZ8W5Ws2IZaERknIhsFpEtInLSqR9E5DoRKRGRj0Xk+ebq1MFgpaKMHDmSiRMnMnToUILBIEOGDOGWW5xfoxQP9WM0sSAiAeAxYCxQBrwnIkuNMSVRZbKAXwDfMcYcEJFmZ9bTQKNUA3PnzmXu3LnxboZj1qFTzL7KI4AtxpgvAESkABgPlESV+TfgMWPMAQBjzJ7mKtVAo5TPWYPBIafFu4tI9JWLC40xC6OeZwA7op6XASMb1HEugIisxbq33z3GmCYvt06KXKdoXk8g5z7vx11eEUBKMIU6F+9JCQp1Yee7ym39XudGteg9mut0XHZue/NMUX9HZUfIpiZznURkIjDOGJNvP78BGGmMmRFVZhlQC1wHZAJvA4OMMQcbqzfpejSmNswI83fH5TfIxYwxyxyXXy1Xu558zU1eEVi5RcvMGMflr5bV7DftHZfvJpWu63ebG+XmMwLrc3K7H9zu5+QV00OncqBX1PNM+7VoZcC7xphaoFREPgWygPcaq1TPOinlczE+vf0ekCUifUUkBEwCljYo8wpwCYCIdMc6lPqiqUqTrkej1KkoVtfIGGPCIjIDeB1r/OUpY8zHInIvUGSMWWovu1xESoAI8DNjzL6m6tVAo5TPxTrXyRizHFje4LU5Uf83wB32wxENNEr5nEGoTvBpEDTQKOVzfsje1kCjlM9poFFKtYl43gLCCQ00SvlcjFMQPJHYrVNKNUsPnTxy4MABKisr2bJly/FpVlJSUnQiOXVKss46Oc51igtf5jodPnyYd999l/POO49PPvkEYwzGGOrq6lzP6+R17pLbvCKAQFCIuMhdCgYh7CKNx239XudGWe/xPtdp2IWDE3W6lVb9QvbJPd3MKRrnqOx0eV7ndXIqFAoRDAY544wz+Pzzz09YZmrDrucHutwUOi7/hoznJvMHx+UXya2u5lwCa94lt7lLNac5rz90yLiu3+28UW4+I7A+J7f7we1+djrdysGDB8nPz2fTpk2ICE899RSjRo1yvK540EMnpXxm5syZjBs3jiVLllBTU0NFRUW8m9QkHaNRymcOHTrE22+/zaJFiwCr9xwKJfb4hx8CjWZvKxWltLSUHj16cPPNNzNkyBDy8/M5duyYo/fOmTOH+fPnH39+55138sgjj3jU0q/VpyA4ecSLBhqlooTDYTZu3Mitt95KcXExHTt2ZN68eY7eO23aNJ555hkA6urqKCgoYOrUqV42F9BZEJTynczMTDIzMxk50rp75cSJE9m4caOj9/bp04fTTz+d4uJi3njjDYYMGcLpp5/uZXOPS/RAo2M0SkU588wz6dWrF5s3b+a8885j9erVDBgwwPH78/PzWbRoEV9++SXTpk3zsKVfi+UsCF7RQKNUA7///e+ZMmUKNTU19OvXj6efftrxe6+55hrmzJlDbW0tzz/f7HRHMaEpCEr50ODBgx1fc9NQKBTi0ksvpWvXrgQCbdfLSPSzThpolIqhuro61q9fz8svv9xm63Q53Upc6GCwUjFSUlJC//79GTNmDFlZWW223voxGiePePFlrlMs53XyPtfJ3ZxL4D53KQi4mbHIbf1ut6Etcp1aUr4uSed1OiM301xfNKP5gsB8+YXmOsWCzuvUPJ3XKfnoGI1SylN+SEHQQKOUz+l1NEopz1lnnXS6FaWUh/TQSSnVJjTQKKU8pWM0SinPaa6TUspzfkhB0ECjlM/poVMbqaurIxwOE4m4mH5DqSSS6IdOvsx1OnLkCOvWraN79+6Ul5eTkpJCIBAgGAxy+VVXJliuk87r5Ow9iTWvUyQSITc3l4yMDJYtc54a0UKtynXqkptlRhTNd1R2tVytuU5OHTlyhLq6OtLT0zl06NAJy0xtmMFmneO6PpBRrucHutL82XH5FTKB6WaB4/IAT8oM17lFbuddclu/m214Uma4+ozA+pzc7ge3+9nNPWYeeeQRsrOzOXz4sOP3xIsfrqPx5W0iunXrRrt27ejWrVu8m6KSUFlZGa+99hr5+fnxbopjiX6bCF/2aJTy0u23384DDzzAkSNH4t0UR+pISfgUBF/2aJTyyrJly+jZsyfDhg2Ld1NcieUsCCIyTkQ2i8gWEZndRLkJImJEpNkxH+3RKBVl7dq1LF26lOXLl1NVVcXhw4eZOnUqzz33XLyb1qhYjtGISAB4DBgLlAHvichSY0xJg3KdgZnAu07q1R6NUlHuv/9+ysrK2Lp1KwUFBYwePTqhgwxYp4ZjOEYzAthijPnCGFMDFADjT1Lu18B/AVVOKtVAo5TvWSkITh5AdxEpino0PF2ZAeyIel5mv/b12kSGAr2MMa85baEeOinViEsuuYRLLrkk3s1olstDp72tuY5GRFKA3wE3uXmfBhqlfM4gVMcu16kc6BX1PNN+rV5nIAd4S0QAzgSWikieMabRC5U00CjlczHO3n4PyBKRvlgBZhJw/fF1GXMI6F7/XETeAmY1FWRAA41SSSFWZ52MMWERmQG8DgSAp4wxH4vIvUCRMWZpS+rVQKOUz8U6BcEYsxxY3uC1OY2UvcRJnb5MqozlBHJuk/O8TsIE90mMbid48zpJ0nWCZEve04KkSld/F22rVUmVwdwLTdcNy5svCOwLZGpSZSyY2jDZZqPj8p/I0IRKwgQrwdDtJHU3mT84Lr9IbnVdv9tEUjefEbQsSdLtfk5Wpk6orkrsFISkCzRKnWqMESLhxM7e1kCjlN8ZNNAopbxljBCu1UCjlPKUUBdJ7K9yYrdOKdU8A+ihk1LKU3UCVYn9VU7s1imlnEnYS4QsSRVodLoVdUqybkiT0JLifjTGGKqqqqiqcnQPHnWKys3NZdy4cfFuRuzVBxonjzjxdY/GGENtbS3V1dW0a9eOdu0S++pIFV9OplvZsWMHP/nJT9i9ezciwi233MLMmTPboHWtYIDaeDeiab7NdVq7di3t2rXj8OHDpKWlYd8bw/NcJ69zo6z3uM0t8rq82/wu99vs9X5wmuu0a9cudu3axdChQzly5AjDhg3jlVdeYcCAAc7X5V6rcp3k/FzDQodzVl0smuvk1J49e6isrCQnJ4eSkhPumYypDXOO2eS4rs8lh3PNh47LfyoXepobBVbezxjjfHbE1XI1l5tCx+XfkPGu6/dycjewPie3+8HtfnYiPT2d9PR0ADp37kx2djbl5eVeB5rW8cEYjS8DTbdu3ejYsSNdu3aNd1NUEtu6dSvFxcWMHDky3k1pmgYabwSDvmy28pGjR48yYcIE5s+fT5cuXeLdnKZpoFHKf2pra5kwYQJTpkzh2muvjXdzmqeBRil/McYwffp0srOzueOOO+LdHOcSPNAkxXU0SsXK2rVrefbZZ3nzzTcZPHgwgwcPZvlyZ3evi5s6rGncnDziRHs0SkX57ne/SzOXfCQePXRSSnlOA41SynMaaJRSbUIDjVLKU9qjaXuSGnR8uTkAwQCfyoWuyn8go1yVXyejnZfHyhVaLVe7Kv+GjPe0flfb4PYzst/jdj+43c9Jqw6ojHcjmpZ0gcbUhjnbfOK4/DbJdp0z4+W8UWDl/Ywwf3dcfoNc7DoXyW39Xs65BNbn5HY/uN3PScsACX4rpqQLNEqdkvTQSSnlKR2jUUp5TgONUspz9SkICUwDjVLJQHs0SilP6aGTUspzPrg5edLcJsIYQ0VFRbyboRJYUk+3EnH4iJOkCDT1QUZv8amaUlRUxMqVK5stt3LlSs477zz69+/PvHnz2qBlrRTjeZ1EZJyIbBaRLSIy+yTL7xCREhH5p4isFpGzm6vT94GmPsiEQiFSU1Pj3Rzlc5FIhNtuu40VK1ZQUlLCCy+88I2ZNhKOwUpBcPJohogEgMeAK4EBwGQRaTgFRDGQa4y5AFgCPNBcvb7uAtTU1FBRUUG7du0IBoMYY+jUtaury83d5kZJapBPZKir8m7zfiQ1yAa52FV5N7lILanfzTa4/Yzq3+N2P7jdz0eOHKFz585NltuwYQP9+/enX79+AEyaNInCwsLEn24ldodFI4AtxpgvAESkABgPHI+2xpj/iSq/HpjaXKW+DTR1dXUUFRWdEGQCgQBLCgqora2lffv2xyeVi6Xa2loikQhpaWkxrxus4CkinvXOqqqqSE1NJRDwJsmwoqKCtLQ0UlJi31mu7722b9/eVf2zZ89m37599O/fn4yMDHr27NnoIVR5eTm9evU6/jwzM5N333231W33lLuzTt1FJHq2uYXGmIVRzzOAHVHPy4Cm5puZDqxobqW+PHSqrKykoqKC888//4QgU11d7WmQMcZQU1Pj6dS74XDYsyAAICKe3qoyFApRU1PjSd0iQlpaGpWVla62Yd68efzpT38iPT2d8vJy9uzZk1yDwu7GaPYaY3KjHgtPWqcDIjIVyAUebK6sLwPNkSNHaN++Pd26dTshyFRXV3sWZACqq6sJhUKe1W+MwRjjSW+gnohQV+d8Oly3gsEgkUjEs2AWCAQIhUJUVbm/FPZ3v/tds8EmIyODHTu+/kEvKysjIyOjVW32XP3pbSeP5pUDvaKeZ9qvnUBELgPuBPKMMdXNVerLQNOzZ08CgQCRSOSEINOhQwfPgkAkEqGurs7TAef67fFSSkqK5zffTk1NpbbWuws7UlNTEZEW9ZyaCzbDhw/ns88+o7S0lJqaGgoKCsjLy4tV070Tu9Pb7wFZItJXRELAJGBpdAERGQL8ESvI7HFSqS/HaOp/+cPhMP369ePTTz/loosu8uyQpq6ujuLiYgYNGkSHDh08WQdAaWkpHTt2pGfPnp6t4+DBg+zZs4dzzz3Xs3WEw2GKi4sZNmyYp4H/ww8/5JxzzuG0005z9d63336bF154gRdffJGXX375hGXBYJAFCxZwxRVXEIlEmDZtGgMHDoxl02MvhrlOxpiwiMwAXgcCwFPGmI9F5F6gyBizFOtQqRPwsr1/txtjmozG0syvW0LOO1FVVcX69eupq6sjHA6TkpLi6eFGfVDz+vR5/fiMV19OsLYlEol4fs1RW2xLXV0ddXV1LdqW2bNn89VXX2GMoXv37nTv3t3RNTYeadWHJB1zDQOKmi8IUCTvG2NyW7O+lvBloFEqybQu0HTINZzvMNAUxyfQ+PLQSSnVgCZVKqU8pdnbSinP6Y2vlFKe0x6NUqpNaKBRSnnKBze+0kCjlN/5YAI5X6YgKBVL+/fvZ+zYsWRlZTF27FgOHDhw0nKLFy8mKyuLrKwsFi9e/I3leXl55OS4mKY3VmJ84ysvaKBRp7x58+YxZswYPvvsM8aMGXPSu+rt37+fuXPn8u6777Jhwwbmzp17QkD6y1/+QqdOndqy2V+rn3s7Bje+8ooGGnXKKywsJC8vj7Fjx7Jw4UIWLFjwjV7N66+/ztixY3n11VcZMWIE+/bt45e//CUAe/bsIT8/n/Xr17NlyxZmz/7G3S+9p/cMViqx7d69m1/96ld88MEHpKSkUF1d/Y1ezbZt21izZg35+fl06dKFKVOm8NJLL3HgwAHGjx9/PKerZ8+erF27lhUrmr0XVGwZh4840UCjTgmXXXYZOTk533gUFhZijKGwsJBXX32VkpISjDG8+OKLJ7x/w4YNVFdXM23aNH72s5+xZs0asrKymDt3LiUlJezcuZPFixfz5ZdfMnjwYMrKyuK0pYlJkyrVKa93797s3LmTZcuWcdttt7F161ZEhHD469HTCy64gNraWoLBIJWVlXzxxRdcccUVbN++ndLS0uP3KwqHwwSDQTZv3nz8vsMOtC6pUnINOEyqJD5JldqjUaeMxno1Z555JsYYbrvtNiZMmMCQIUOIRCL07t37+CFUbW0tW7duJTU1lX/84x8YY1i9ejXbt2/n1ltvpbS0lE8++YRAIEAgEOCCCy7gt7/9bZy3OHHodTTqlLFq1aqTvv7000+Tn5/Pzp07ef/99yktLSUUCnHVVVfx4IMPkpeXRyAQoGvXrmzevJlBgwaRlpaGMYY+ffqwdu1aHnroIa677jpSUlIIh8N8//vfb8Mtqz/tlLi0R6NOeeeffz6nn346OTk5/OY3v0FE6NOnD8FgkFmzZlFYWEhGRgbt2rWjW7du5OTkkJqaSiQS4bLLLuOjjz7izjvvJBKJcNZZZ9G9e3cGDRrUhlsQ25sGe0EDjTrlDR8+HGMM27dv55prruHQoUP8+Mc/BqxDpmeffZa8vDwOHz7MzJkzWbNmDeFwmG9961tMmDCB2tpa7rvvPtavX8+2bdvo0qULGzdubMMtSPwr9jTQqFNeMBhk9uzZHD16FGMMF154Ie3bt6ekpIQDBw4wevRopk+fTigUYv78+YRCIVauXEldXR3f+c536NSpExkZGXz55ZdMnTqVTz/9lKFD3U2g1zqJ36PRMRqlgJkzZ/LYY4/x4IMP8vjjj1NQUMDzzz/P0qXWBABpaWncfffdfPTRR5SWlvLWW28xevRoIpEIIkK3bt3o0qULa9asoU+fPhw8eJCUlBTS0tKYMWOGx61P/KxKDTRK8fXsBzNnzmTr1q389Kc/JSsri/nz53P33XcDMH36dG644QY2btxIcXEx69evZ8mSJaSnp/P555+TlZVFKBQC4Ec/+hFnnHFGGwQZ+Hry7cSl19Eo1cDy5cu5/fbbj0+3cueddzJnzhxyc3PJy8ujqqqKG264geLiYrp160ZBQcE3rpm555576NSpE7NmzXKyylZeR5Nj4OXmCwIwQGdBUOoU1cpAM9DACw5LX6izICilWiLx7+WpgUYp30v8wWA9va1UG3j88ccZPHgwgwcPpm/fvlx66aUxrD3xr6PRMRql2lBtbS2jR4/m5z//OT/4wQ/qX27lGM25Bv7bYemxOkajVLKbOXMmo0ePjg4yMZD4h04aaJRqI4sWLWLbtm0sWLDAg9p1MFipU97777/Pb3/7W9asWUNKSqyHRrVHo5QCFixYwP79+48PAufm5vLEE0/EqHYNNEoprHveeEevo1FKeS7xb3ylgUYp39NDJ6WU5xL/0EmvDFbK92J74ysRGScim0Vki4h8YzY8EWknIi/ay98VkT7N1amBRinfi10KgogEgMeAK4EBwGQRGdCg2HTggDGmP/Aw8F/N1auBRinfi+nk2yOALcaYL4wxNUABML5BmfHAYvv/S4AxUj9VZyOaG6NpVQ6GUqot7Hod7unusHCaiETPNrfQGLMw6nkGsCPqeRkwskEdx8sYY8Iicgg4Hdjb2Ep1MFgpnzPGjIt3G5qjh05KqWjlQK+o55n2ayctIyJB4DRgX1OVaqBRSkV7D8gSkb4iEgImAUsblFkK3Gj/fyLwpmnmfjN66KSUOs4ec5kBvA4EgKeMMR+LyL1AkTFmKfAk8KyIbAH2YwWjJjV34yullGo1PXRSSnlOA41SynMaaJRSntNAo5TynAYapZTnNNAopTyngUYp5bn/D3G/GbNuR1qfAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# define iteration parameters\n",
    "newton_tol = 1e-6\n",
    "t = .0\n",
    "T = 1\n",
    "k = 0\n",
    "times = [t]\n",
    "\n",
    "# Time loop\n",
    "while t < T:\n",
    "    # Increment time\n",
    "    t += dt\n",
    "    k += 1\n",
    "    times.append(t)\n",
    "    p0 = p.val\n",
    "    print('Solving time step: ', k)\n",
    "    # solve newton iteration\n",
    "    err = np.inf\n",
    "    while err > newton_tol:\n",
    "        eq = f(p, p0)   \n",
    "        p = p - sps.linalg.spsolve(eq.jac, eq.val)\n",
    "        err = np.sqrt(np.sum(eq.val**2))\n",
    "    pp.plot_grid(g, p.val,color_map = [0, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
